From 51cbf02119bf64027e774b5e22573791db67ec90 Mon Sep 17 00:00:00 2001 From: Andrew Howard Date: Sun, 10 Dec 2023 21:19:02 -0800 Subject: [PATCH 001/153] Added recipe for QC development checklist --- docs/source/develop/start.rst | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/docs/source/develop/start.rst b/docs/source/develop/start.rst index 95eeb6277..6a59a70db 100644 --- a/docs/source/develop/start.rst +++ b/docs/source/develop/start.rst @@ -21,7 +21,7 @@ Here are the steps to adding a new quality control test. #. Start a Git branch for your feature. #. Write a method for your QC check in `KPF-modules/quality_control/src/quality_control.py `_ based on code from your Jupyter notebook. The method should return a boolean (``QC_pass``) that is True if the input KPF object passed the QC check and False otherwise. One method to model yours on is ``L0_data_products_check()``. Your method should be in the appropriate class for the data level of your QC check. For example, for a QC check to an L0 object, put the method in the ``QCL0`` class in ``quality_control.py``. #. Add information about your QC to the QCDefinitions class in ``quality_control.py``. You can model your dictionary entries on the ones for ``name4 = 'L0_data_products_check'``. -#. Check that your QC works as expected. See `this Jupyter notebook `_ for an example. +#. Check that your QC works as expected. See `this Jupyter notebook `_ for an example. You can also modify the config file specified in this command and check the result: ``kpf -c configs/quality_control_example.cfg -r recipes/quality_control_example.recipe``. #. Commit the changes to your Git branch and submit a pull request. Developing Quicklook Plots From bfe70cf266bd2bdc8fa13152e3b40f0e18086b10 Mon Sep 17 00:00:00 2001 From: Andrew Howard Date: Sun, 10 Dec 2023 21:27:39 -0800 Subject: [PATCH 002/153] Added instruction re: keyword documentation --- docs/source/develop/start.rst | 1 + 1 file changed, 1 insertion(+) diff --git a/docs/source/develop/start.rst b/docs/source/develop/start.rst index 6a59a70db..2648d1fb3 100644 --- a/docs/source/develop/start.rst +++ b/docs/source/develop/start.rst @@ -23,6 +23,7 @@ Here are the steps to adding a new quality control test. #. Add information about your QC to the QCDefinitions class in ``quality_control.py``. You can model your dictionary entries on the ones for ``name4 = 'L0_data_products_check'``. #. Check that your QC works as expected. See `this Jupyter notebook `_ for an example. You can also modify the config file specified in this command and check the result: ``kpf -c configs/quality_control_example.cfg -r recipes/quality_control_example.recipe``. #. Commit the changes to your Git branch and submit a pull request. +#. Document the new QC-related FITS keywords int the appropriate section of 'KPF Data Format' in Readthedocs. Developing Quicklook Plots ^^^^^^^^^^^^^^^^^^^^^^^^^^ From afc116e269a40b6944607955f2ca49a45bc587ea Mon Sep 17 00:00:00 2001 From: Andrew Howard Date: Sun, 10 Dec 2023 21:33:12 -0800 Subject: [PATCH 003/153] added exposure_meter_not_saturated_check QC method --- configs/quality_control_example.cfg | 4 +- .../quality_control/src/quality_control.py | 77 +++++++++++++++++++ modules/quicklook/src/analyze_em.py | 8 +- 3 files changed, 83 insertions(+), 6 deletions(-) diff --git a/configs/quality_control_example.cfg b/configs/quality_control_example.cfg index bf870720a..8c3ce38a4 100644 --- a/configs/quality_control_example.cfg +++ b/configs/quality_control_example.cfg @@ -20,8 +20,8 @@ data_type = KPF #output_fits_filename = /testdata/L1/20231030/KP.20231030.00336.79_L1.fits data_level_str = L0 -input_fits_filename = /data/L0/20230829/KP.20230829.76026.81.fits -output_fits_filename = /testdata/L0/20230829/KP.20230829.76026.81.fits +input_fits_filename = /data/L0/20230701/KP.20230701.79161.47.fits +output_fits_filename = /testdata/L0/20230701/KP.20230701.79161.47.fits [MODULE_CONFIGS] quality_control = modules/quality_control/configs/default.cfg diff --git a/modules/quality_control/src/quality_control.py b/modules/quality_control/src/quality_control.py index 7e3630478..49fc60e9a 100644 --- a/modules/quality_control/src/quality_control.py +++ b/modules/quality_control/src/quality_control.py @@ -148,6 +148,16 @@ def __init__(self): self.db_columns[name4] = None self.methods[name4] = ["add_qc_keyword_to_header"] + name5 = 'exposure_meter_not_saturated_check' + self.names.append(name5) + self.kpf_data_levels[name5] = ['L0'] + self.descriptions[name5] = 'Check if 2+ reduced EM pixels are within 90% of saturation in EM-SCI or EM-SKY.' + self.data_types[name5] = 'int' + self.fits_keywords[name5] = 'EMSAT' + self.fits_comments[name5] = 'QC: EM not saturated check' + self.db_columns[name5] = None + self.methods[name5] = ["add_qc_keyword_to_header"] + # Integrity checks. if len(self.names) != len(self.kpf_data_levels): raise ValueError("Length of kpf_data_levels list does not equal number of entries in descriptions dictionary.") @@ -441,6 +451,73 @@ def L0_header_keywords_present_check(self, essential_keywords=['auto'], debug=Fa return QC_pass + def exposure_meter_not_saturated_check(self, debug=False): + """ + This Quality Control function checks if 2 or more reduced pixels in an exposure + meter spectrum is within 90% of saturated. The check is applied to the EM-SCI + and EM-SKY fibers and returns False if saturation is detected in either. + Note that this check only works for L0 files with the EXPMETER_SCI and + EXPMETER_SKY extensions present. + + Args: + L0 - an L0 object + fiber ('SCI' [default value] or 'SKY) - the EM fiber output to be tested + debug - an optional flag. If True, missing data products are noted. + + Returns: + QC_pass - a boolean signifying that the QC passed (True) for failed (False) + """ + + saturation_level = 1.93e6 # saturation level in reduced EM spectra (in data frame) + saturation_fraction = 0.9 + + # Read and condition the table of Exposure Meter Data + L0 = self.kpf_object + if hasattr(L0, 'EXPMETER_SCI') and hasattr(L0, 'EXPMETER_SKY'): + if (L0['EXPMETER_SCI'].size > 1) and (L0['EXPMETER_SKY'].size > 1): + pass + else: + return False + else: + return True # pass test if no exposure meter data present + EM_sat_SCI = L0['EXPMETER_SCI'].copy() + EM_sat_SKY = L0['EXPMETER_SKY'].copy() + columns_to_drop_SCI = [col for col in EM_sat_SCI.columns if col.startswith('Date')] + columns_to_drop_SKY = [col for col in EM_sat_SKY.columns if col.startswith('Date')] + EM_sat_SCI.drop(columns_to_drop_SCI, axis=1, inplace=True) + EM_sat_SKY.drop(columns_to_drop_SKY, axis=1, inplace=True) + if len(EM_sat_SCI) >= 3: # drop first and last rows if nrows >= 3 + EM_sat_SCI = EM_sat_SCI.iloc[1:-1] + EM_sat_SKY = EM_sat_SKY.iloc[1:-1] + + # Determine the saturation fraction + for col in EM_sat_SCI.columns: + try: # only apply to columns with wavelengths as headers + float_col_title = float(col) + EM_sat_SCI[col] = EM_sat_SCI[col] / saturation_level + except ValueError: + pass + for col in EM_sat_SKY.columns: + try: + float_col_title = float(col) + EM_sat_SKY[col] = EM_sat_SKY[col] / saturation_level + except ValueError: + pass + + saturated_elements_SCI = (EM_sat_SCI > saturation_fraction).sum().sum() + saturated_elements_SKY = (EM_sat_SKY > saturation_fraction).sum().sum() + total_elements = EM_sat_SCI.shape[0] * EM_sat_SCI.shape[1] + saturated_fraction_threshold = 1.5 / EM_sat_SCI.shape[1] + + if saturated_elements_SCI / total_elements > saturated_fraction_threshold: + QC_pass = False + elif saturated_elements_SKY / total_elements > saturated_fraction_threshold: + QC_pass = False + else: + QC_pass = True + + return QC_pass + ##################################################################### class QC2D(QC): diff --git a/modules/quicklook/src/analyze_em.py b/modules/quicklook/src/analyze_em.py index 63e14c643..bf900c616 100644 --- a/modules/quicklook/src/analyze_em.py +++ b/modules/quicklook/src/analyze_em.py @@ -200,11 +200,11 @@ def plot_EM_time_series(self, fiber='both', fig_path=None, show_plot=False): plt.plot(self.time_em, self.flux_SKY /self.SKY_SCI_ratio/ self.flux_SCI , marker='o', linewidth=4, color='k', label = r'SKY$_{\mathrm{corrected}}$ / SCI - 445-870 nm') plt.plot(0, 0, marker='o', markersize=0.1, color='white') # force the y-axis to go to zero if fiber == 'both': - plottitle = 'Exposure Meter Time Series (SCI and SKY) - ' + str(self.EM_nexp)+ ' EM exposures, ' + str(self.EM_texp)+ ' sec - ' + str(self.ObsID) + ' - ' + self.name + plottitle = 'Exposure Meter Time Series (SCI and SKY) - ' + str(self.EM_nexp)+ r' EM exposures $\times$ ' + str(self.EM_texp)+ ' sec - ' + str(self.ObsID) + ' - ' + self.name elif fiber == 'sky': - plottitle = 'Exposure Meter Time Series (SKY) - ' + str(self.EM_nexp)+ ' EM exposures, ' + str(self.EM_texp)+ ' sec - ' + str(self.ObsID) + ' - ' + self.name + plottitle = 'Exposure Meter Time Series (SKY) - ' + str(self.EM_nexp)+ r' EM exposures $\times$ ' + str(self.EM_texp)+ ' sec - ' + str(self.ObsID) + ' - ' + self.name elif fiber == 'sci': - plottitle = 'Exposure Meter Time Series (SCI) - ' + str(self.EM_nexp)+ ' EM exposures, ' + str(self.EM_texp)+ ' sec - ' + str(self.ObsID) + ' - ' + self.name + plottitle = 'Exposure Meter Time Series (SCI) - ' + str(self.EM_nexp)+ r' EM exposures $\times$ ' + str(self.EM_texp)+ ' sec - ' + str(self.ObsID) + ' - ' + self.name elif fiber == 'ratio': median_flux_ratio = np.nanmedian(self.flux_SKY/self.SKY_SCI_ratio / self.flux_SCI) avg_flux_ratio = np.nansum(self.flux_SKY/self.SKY_SCI_ratio) / np.nansum(self.flux_SCI) @@ -290,7 +290,7 @@ def plot_EM_spectrum(self, fig_path=None, show_plot=False): lns2 = ax2.plot(self.wav_SKY, self.int_SKY_spec, marker='.', color='brown', label = 'SKY',zorder = 0, alpha = 0.5) ax1.set_ylim(0,np.nanpercentile(self.int_SCI_spec,99.9)*1.1) ax2.set_ylim(0,np.nanpercentile(self.int_SKY_spec,99.9)*1.1) - plt.title('Exposure Meter Spectrum: ' + str(self.EM_nexp)+ ' EM exposures, ' + str(self.EM_texp)+ ' sec - ' + str(self.ObsID) + ' - ' + self.name, fontsize=14) + plt.title('Exposure Meter Spectrum: ' + str(self.EM_nexp)+ r' EM exposures $\times$ ' + str(self.EM_texp)+ ' sec - ' + str(self.ObsID) + ' - ' + self.name, fontsize=14) plt.yticks(fontsize=14, color='brown') ax1.set_xlabel("Wavelength (nm)",fontsize=14) ax1.set_ylabel("Exposure Meter SCI Flux (e-/nm/s)",fontsize=14) From 4db65ac915f7870143a74f7c691b9b795da7e4b2 Mon Sep 17 00:00:00 2001 From: Andrew Howard Date: Sun, 10 Dec 2023 23:05:46 -0800 Subject: [PATCH 004/153] added QC method exposure_meter_flux_not_negative_check() --- configs/quality_control_example.cfg | 4 +- .../quality_control/src/quality_control.py | 101 +++++++++++++++--- 2 files changed, 86 insertions(+), 19 deletions(-) diff --git a/configs/quality_control_example.cfg b/configs/quality_control_example.cfg index 8c3ce38a4..19e68f79f 100644 --- a/configs/quality_control_example.cfg +++ b/configs/quality_control_example.cfg @@ -20,8 +20,8 @@ data_type = KPF #output_fits_filename = /testdata/L1/20231030/KP.20231030.00336.79_L1.fits data_level_str = L0 -input_fits_filename = /data/L0/20230701/KP.20230701.79161.47.fits -output_fits_filename = /testdata/L0/20230701/KP.20230701.79161.47.fits +input_fits_filename = /data/L0/20231204/KP.20231204.43445.39.fits +output_fits_filename = /testdata/L0/20231204/KP.20231204.43445.39.fits [MODULE_CONFIGS] quality_control = modules/quality_control/configs/default.cfg diff --git a/modules/quality_control/src/quality_control.py b/modules/quality_control/src/quality_control.py index 49fc60e9a..8f67edca1 100644 --- a/modules/quality_control/src/quality_control.py +++ b/modules/quality_control/src/quality_control.py @@ -2,6 +2,7 @@ import numpy as np import numpy.ma as ma import pandas as pd +from scipy.ndimage import convolve1d from modules.Utils.kpf_parse import get_data_products_L0 """ @@ -81,7 +82,7 @@ class QCDefinitions: Must be a short string for brevity (say, under 35 characters), following the FITS standard. db_columns (dictionary of strings): Each entry specifies either database_table.column if applicable, or None if not. - methods (dictionary of lists): Each entry specifies a list of methods that apply to the metric. +# methods (dictionary of lists): Each entry specifies a list of methods that apply to the metric. """ def __init__(self): @@ -93,7 +94,7 @@ def __init__(self): self.fits_keywords = {} self.fits_comments = {} self.db_columns = {} - self.methods = {} + #self.methods = {} # Define QC metrics @@ -105,7 +106,7 @@ def __init__(self): self.fits_keywords[name0] = 'JBTRED1' self.fits_comments[name0] = 'QC: J-B test for RED AMP-1 detector' self.db_columns[name0] = None - self.methods[name0] = ["add_qc_keyword_to_header"] + #self.methods[name0] = ["add_qc_keyword_to_header"] name1 = 'not_junk_check' self.names.append(name1) @@ -115,7 +116,7 @@ def __init__(self): self.fits_keywords[name1] = 'NOTJUNK' self.fits_comments[name1] = 'QC: Not in list of junk files check' self.db_columns[name1] = None - self.methods[name1] = ["add_qc_keyword_to_header"] + #self.methods[name1] = ["add_qc_keyword_to_header"] name2 = 'monotonic_wavelength_solution_check' self.names.append(name2) @@ -126,7 +127,7 @@ def __init__(self): self.fits_comments[name2] = 'QC: Monotonic wavelength-solution check' self.db_columns[name2] = None #self.methods[name2] = ["add_qc_keyword_to_header","monotonic_wavelength_solution_check","add_qc_keyword_to_header_for_monotonic_wls"] - self.methods[name2] = ["add_qc_keyword_to_header"] + #self.methods[name2] = ["add_qc_keyword_to_header"] name3 = 'L0_data_products_check' self.names.append(name3) @@ -136,7 +137,7 @@ def __init__(self): self.fits_keywords[name3] = 'DATAPRL0' self.fits_comments[name3] = 'QC: L0 data present check' self.db_columns[name3] = None - self.methods[name3] = ["add_qc_keyword_to_header"] + #self.methods[name3] = ["add_qc_keyword_to_header"] name4 = 'L0_header_keywords_present_check' self.names.append(name4) @@ -146,7 +147,7 @@ def __init__(self): self.fits_keywords[name4] = 'KWRDPRL0' self.fits_comments[name4] = 'QC: L0 keywords present check' self.db_columns[name4] = None - self.methods[name4] = ["add_qc_keyword_to_header"] + #self.methods[name4] = ["add_qc_keyword_to_header"] name5 = 'exposure_meter_not_saturated_check' self.names.append(name5) @@ -156,9 +157,19 @@ def __init__(self): self.fits_keywords[name5] = 'EMSAT' self.fits_comments[name5] = 'QC: EM not saturated check' self.db_columns[name5] = None - self.methods[name5] = ["add_qc_keyword_to_header"] - - # Integrity checks. + #self.methods[name5] = ["add_qc_keyword_to_header"] + + name6 = 'exposure_meter_flux_not_negative_check' + self.names.append(name6) + self.kpf_data_levels[name6] = ['L0'] + self.descriptions[name6] = 'Check for negative flux in the EM-SCI and EM-SKY by looking for 2o0 consecuitive pixels in the summed spectra with negative flux.' + self.data_types[name6] = 'int' + self.fits_keywords[name6] = 'EMNEG' + self.fits_comments[name6] = 'QC: EM not negative flux check' + self.db_columns[name6] = None + #self.methods[name6] = ["add_qc_keyword_to_header"] + + # Integrity checks if len(self.names) != len(self.kpf_data_levels): raise ValueError("Length of kpf_data_levels list does not equal number of entries in descriptions dictionary.") @@ -177,8 +188,8 @@ def __init__(self): if len(self.names) != len(self.db_columns): raise ValueError("Length of db_columns list does not equal number of entries in db_columns dictionary.") - if len(self.names) != len(self.methods): - raise ValueError("Length of methods list does not equal number of entries in methods dictionary.") +# if len(self.names) != len(self.methods): +# raise ValueError("Length of methods list does not equal number of entries in methods dictionary.") keys_list = self.data_types.keys() for key in keys_list: @@ -190,7 +201,7 @@ def __init__(self): def list_qc_metrics(self): - print("name | data_type | keyword | comment | methods | db_column | description |") + print("name | data_type | keyword | comment | db_column | description |") qc_names = self.names @@ -199,11 +210,11 @@ def list_qc_metrics(self): data_type = self.data_types[qc_name] keyword = self.fits_keywords[qc_name] comment = self.fits_comments[qc_name] - methods = self.methods[qc_name] +# methods = self.methods[qc_name] db_column = self.db_columns[qc_name] description = self.descriptions[qc_name] - print(qc_name," | ",data_type," | ",keyword," | ",comment," | ",methods," | ",db_column," | ",description) + print(qc_name," | ",data_type," | ",keyword," | ",comment," | ",db_column," | ",description) ##################################################################### @@ -351,8 +362,8 @@ def L0_data_products_check(self, debug=False): L0 = self.kpf_object - # determine which extensions should be in the L0 file - # first add triggrered cameras (Green, Red, CaHK, ExpMeter) + # Determine which extensions should be in the L0 file. + # First add the triggrered cameras (Green, Red, CaHK, ExpMeter) to list of data products trigtarg = L0.header['PRIMARY']['TRIGTARG'] if len(trigtarg) > 0: data_products = trigtarg.split(',') @@ -518,6 +529,62 @@ def exposure_meter_not_saturated_check(self, debug=False): return QC_pass + + def exposure_meter_flux_not_negative_check(self, debug=False): + """ + This Quality Control function checks if 20 or more consecutive elements of the + exposure meter spectra are negative. Negative flux usually indicates + over-subtraction of bias from the raw EM images. The check is applied to the + EM-SCI and EM-SKY fibers and returns False if negative flux is detected in + either. Note that this check only works for L0 files with the EXPMETER_SCI and + EXPMETER_SKY extensions present. + + Args: + L0 - an L0 object + fiber ('SCI' [default value] or 'SKY) - the EM fiber output to be tested + debug - an optional flag. If True, missing data products are noted. + + Returns: + QC_pass - a boolean signifying that the QC passed (True) for failed (False) + """ + + N_in_a_row = 20 # number of negative flux elements in a row that triggers QC failure + + # Read and condition the table of Exposure Meter Data + L0 = self.kpf_object + if hasattr(L0, 'EXPMETER_SCI') and hasattr(L0, 'EXPMETER_SKY'): + if (L0['EXPMETER_SCI'].size > 1) and (L0['EXPMETER_SKY'].size > 1): + pass + else: + return False + else: + return True # pass test if no exposure meter data present + EM_SCI = L0['EXPMETER_SCI'].copy() + EM_SKY = L0['EXPMETER_SKY'].copy() + columns_to_drop_SCI = [col for col in EM_SCI.columns if col.startswith('Date')] + columns_to_drop_SKY = [col for col in EM_SKY.columns if col.startswith('Date')] + EM_SCI.drop(columns_to_drop_SCI, axis=1, inplace=True) + EM_SKY.drop(columns_to_drop_SKY, axis=1, inplace=True) + counts_SCI = EM_SCI.sum(axis=0).values + counts_SKY = EM_SKY.sum(axis=0).values + + # Determine if the spectra have significant negative flux + negative_mask_SCI = counts_SCI < 0 # spectral elements with negative flux + negative_mask_SKY = counts_SKY < 0 + window = np.ones(N_in_a_row, dtype=int) # window to convolve with spectra + conv_result_SCI = convolve1d(negative_mask_SCI.astype(int), window, mode='constant', cval=0) + conv_result_SKY = convolve1d(negative_mask_SKY.astype(int), window, mode='constant', cval=0) + has_consec_negs_SCI = np.any(conv_result_SCI == N_in_a_row) + has_consec_negs_SKY = np.any(conv_result_SKY == N_in_a_row) + + if has_consec_negs_SCI or has_consec_negs_SKY: + QC_pass = False + else: + QC_pass = True + + return QC_pass + + ##################################################################### class QC2D(QC): From ecef91ab9d552c02686d3e93f0bae564432d8ed9 Mon Sep 17 00:00:00 2001 From: Aaron Householder Date: Wed, 13 Dec 2023 12:50:56 -0500 Subject: [PATCH 005/153] Current ThAr code --- modules/wavelength_cal/configs/KPF.cfg | 2 +- modules/wavelength_cal/src/alg.py | 120 +++++++++++-------------- 2 files changed, 52 insertions(+), 70 deletions(-) diff --git a/modules/wavelength_cal/configs/KPF.cfg b/modules/wavelength_cal/configs/KPF.cfg index 25fa7bc06..52b31dd54 100644 --- a/modules/wavelength_cal/configs/KPF.cfg +++ b/modules/wavelength_cal/configs/KPF.cfg @@ -13,7 +13,7 @@ fit_order = 4 fit_type = legendre n_sections = 20 clip_peaks = 1 -clip_below_median = 1 +clip_below_median = 0 peak_height_threshold = 0.3 sigma_clip = 2.5 fit_iterations = 5 diff --git a/modules/wavelength_cal/src/alg.py b/modules/wavelength_cal/src/alg.py index 7ecaa0fd3..bee8f0134 100644 --- a/modules/wavelength_cal/src/alg.py +++ b/modules/wavelength_cal/src/alg.py @@ -140,8 +140,8 @@ def run_wavelength_cal( # perform wavelength calibration poly_soln, wls_and_pixels, orderlet_dict = self.fit_many_orders( masked_calflux, order_list, rough_wls=rough_wls, - comb_lines_angstrom=lfc_allowed_wls, - expected_peak_locs=peak_wavelengths_ang, + comb_lines_angstrom=lfc_allowed_wls, + expected_peak_locs=peak_wavelengths_ang, peak_wavelengths_ang=peak_wavelengths_ang, print_update=True, plt_path=self.save_diagnostics_dir ) @@ -172,7 +172,6 @@ def run_wavelength_cal( '{}/all_wls.png'.format(self.save_diagnostics_dir), dpi=250 ) - plt.close() if self.quicklook == True: #TODO @@ -184,8 +183,8 @@ def run_wavelength_cal( poly_soln, wls_and_pixels, orderlet_dict = self.fit_many_orders( masked_calflux, order_list, rough_wls=rough_wls, - comb_lines_angstrom=lfc_allowed_wls, - expected_peak_locs=peak_wavelengths_ang, + comb_lines_angstrom=lfc_allowed_wls, + expected_peak_locs=peak_wavelengths_ang, peak_wavelengths_ang=peak_wavelengths_ang, print_update=True, plt_path=self.save_diagnostics_dir ###CHECK THIS TODO ) @@ -234,7 +233,7 @@ def find_etalon_peaks(self,flux,wave,etalon_mask): def fit_many_orders( self, cal_flux, order_list, rough_wls=None, comb_lines_angstrom=None, - expected_peak_locs=None, plt_path=None, print_update=False): + expected_peak_locs=None,peak_wavelengths_ang=None, plt_path=None, print_update=False): """ Iteratively performs wavelength calibration for all orders. Args: @@ -453,7 +452,7 @@ def fit_many_orders( # calculate the wavelength solution for the order polynomial_wls, leg_out = self.fit_polynomial( - wls, n_pixels, fitted_peak_pixels, peak_heights=peak_heights, + wls, rough_wls_order, peak_wavelengths_ang, order_list, n_pixels, fitted_peak_pixels, peak_heights=peak_heights, plot_path=order_plt_path, fit_iterations=self.fit_iterations, sigma_clip=self.sigma_clip) poly_soln_final_array[order_num,:] = polynomial_wls @@ -662,8 +661,8 @@ def find_peaks(self, order_flux, peak_height_threshold=1.5): detected_peaks, properties = signal.find_peaks(c, distance=distance, height=height) peak_heights = np.array(properties['peak_heights']) - # Only consider peaks with height greater than 5 - valid_peak_indices = np.where(peak_heights > 5)[0] + # Only consider peaks with height greater than 500 + valid_peak_indices = np.where(peak_heights > 500)[0] detected_peaks = detected_peaks[valid_peak_indices] peak_heights = peak_heights[valid_peak_indices] @@ -892,35 +891,6 @@ def clip_peaks( return good_peak_idx - def estimate_fwhm(self, flux, peak_pixel): - """ - Estimate the FWHM of a line at a given pixel position. - - Args: - flux (np.array): The flux values of the spectrum. - peak_pixel (int): The pixel position of the peak of the line. - - Returns: - float: An estimated FWHM value. - """ - half_max = flux[peak_pixel] / 2 - - # Find the points at half maximum on both sides of the peak - left_idx = right_idx = peak_pixel - while left_idx > 0 and flux[left_idx] > half_max: - left_idx -= 1 - while right_idx < len(flux) - 1 and flux[right_idx] > half_max: - right_idx += 1 - - # Interpolate to find more accurate positions where the flux crosses the half maximum - if left_idx > 0: - left_idx += (half_max - flux[left_idx]) / (flux[left_idx + 1] - flux[left_idx]) - if right_idx < len(flux) - 1: - right_idx -= (half_max - flux[right_idx]) / (flux[right_idx - 1] - flux[right_idx]) - - fwhm = right_idx - left_idx - return fwhm - def line_match(self, flux, linelist, line_pixels_expected, plot_toggle, savefig, gaussian_fit_width=10): """ Given a linelist of known wavelengths of peaks and expected pixel locations @@ -1401,6 +1371,11 @@ def fit_gaussian_integral(self, x, y): if self.cal_type == 'ThAr': # Quality Checks for Gaussian Fits + + if max(y) == 0: + print('Amplitude is 0') + return(None, line_dict) + chi_squared_threshold = int(self.chi_2_threshold) # Calculate chi^2 @@ -1431,15 +1406,15 @@ def fit_gaussian_integral(self, x, y): print("Chi squared exceeded the threshold for this line. Line skipped") return None, line_dict - #Check if the Gaussian amplitude is positive and the peak is higher than the wings - if popt[0] <= 0 or popt[0] <= popt[3]: - line_dict['quality'] = 'bad_amplitude' # Mark the fit as bad due to negative amplitude or U shaped gaussian - print('Negative amplitude detected') + # Check if the Gaussian amplitude is positive, the peak is higher than the wings, or the peak is too high + if popt[0] <= 0 or popt[0] <= popt[3] or popt[0] >= 500*max(y): + line_dict['quality'] = 'bad_amplitude' # Mark the fit as bad due to bad amplitude or U shaped gaussian + print('Bad amplitude detected') return None, line_dict return (popt, line_dict) - def fit_polynomial(self, wls, n_pixels, fitted_peak_pixels, fit_iterations=5, sigma_clip=2.1, peak_heights=None, plot_path=None): + def fit_polynomial(self, wls, rough_wls_order, peak_wavelengths_ang, order_list, n_pixels, fitted_peak_pixels, fit_iterations=5, sigma_clip=2.1, peak_heights=None, plot_path=None): """ Given precise wavelengths of detected LFC order_flux lines, fits a polynomial wavelength solution. @@ -1476,15 +1451,22 @@ def fit_polynomial(self, wls, n_pixels, fitted_peak_pixels, fit_iterations=5, si unclipped_idx = np.intersect1d(unclipped_idx, unique_idx[count < 2]) sorted_idx = np.argsort(fitted_peak_pixels[unclipped_idx]) - x, y, w = fitted_peak_pixels[unclipped_idx][sorted_idx], wls[unclipped_idx][sorted_idx], weights[unclipped_idx][sorted_idx] + x, y, w = fitted_peak_pixels[unclipped_idx][sorted_idx], wls[unclipped_idx][sorted_idx], weights[unclipped_idx][sorted_idx] + for i in range(fit_iterations): if self.fit_type.lower() == 'legendre': - leg_out = Legendre.fit(x, y, self.fit_order, w=w) - our_wavelength_solution_for_order = leg_out(np.arange(n_pixels)) + if self.cal_type == 'ThAr': + leg_out = Legendre.fit(x, y, 9, w=w) + our_wavelength_solution_for_order = rough_wls_order # + 1D polynomial, to be added + if self.cal_type == 'lfc': + leg_out = Legendre.fit(x, y, 9, w=w) + our_wavelength_solution_for_order = leg_out(np.arange(n_pixels)) + if self.fit_type == 'spline': leg_out = UnivariateSpline(x, y, w, k=5) our_wavelength_solution_for_order = leg_out(np.arange(n_pixels)) + res = y - leg_out(x) good = np.where(np.abs(res) <= sigma_clip*np.std(res)) @@ -1492,39 +1474,39 @@ def fit_polynomial(self, wls, n_pixels, fitted_peak_pixels, fit_iterations=5, si y = y[good] w = w[good] res = res[good] - # plt.plot(x, res, 'k.') # plt.axhline(0, color='b', lw=2) # plt.xlabel('Pixel') # plt.ylabel('Fit residuals [$\AA$]') # plt.tight_layout() - # plt.savefig('polyfits/polyfit_{}.png'.format(wls[0])) + # plt.savefig('polyfit//polyfit_{}.png'.format(wls[0])) # plt.close() - if plot_path is not None: - sorted_idx = np.argsort(fitted_peak_pixels[unclipped_idx]) - s = InterpolatedUnivariateSpline(fitted_peak_pixels[unclipped_idx][sorted_idx], wls[unclipped_idx][sorted_idx]) - interpolated_ground_truth = s(np.arange(n_pixels)) + #if plot_path is not None: + + # sorted_idx = np.argsort(fitted_peak_pixels[unclipped_idx]) + # s = InterpolatedUnivariateSpline(fitted_peak_pixels[unclipped_idx][sorted_idx], wls[unclipped_idx][sorted_idx]) + # interpolated_ground_truth = s(np.arange(n_pixels)) - approx_dispersion = (our_wavelength_solution_for_order[2000] - our_wavelength_solution_for_order[2100])/100 + # approx_dispersion = (our_wavelength_solution_for_order[2000] - our_wavelength_solution_for_order[2100])/100 # plot ground truth wls vs our wls - fig, ax1 = plt.subplots(tight_layout=True, figsize=(8, 4)) - ax1.plot( - np.arange(n_pixels), - interpolated_ground_truth - our_wavelength_solution_for_order, - color='k' - ) - ax1.set_xlabel('Pixel') - ax1.set_ylabel(r'Wavelength Difference ($\AA$)') - ax2 = ax1.twinx() - warnings.filterwarnings("ignore", "FixedFormatter should only be used together with FixedLocator") - ax2.set_ylabel("Difference (pixels) \nusing dispersion " + r'$\approx$' + '{0:.2}'.format(approx_dispersion) + r' $\AA$/pixel') - ax2.set_ylim(ax1.get_ylim()) - ax1_ticks = ax1.get_yticks() - ax2.set_yticklabels([str(round(tick / approx_dispersion, 2)) for tick in ax1_ticks]) - plt.savefig('{}/interp_vs_our_wls.png'.format(plot_path)) - plt.close() + # fig, ax1 = plt.subplots(tight_layout=True, figsize=(8, 4)) + # ax1.plot( + # np.arange(n_pixels), + # interpolated_ground_truth - our_wavelength_solution_for_order, + # color='k' + # ) + # ax1.set_xlabel('Pixel') + # ax1.set_ylabel(r'Wavelength Difference ($\AA$)') + # ax2 = ax1.twinx() + # warnings.filterwarnings("ignore", "FixedFormatter should only be used together with FixedLocator") + # ax2.set_ylabel("Difference (pixels) \nusing dispersion " + r'$\approx$' + '{0:.2}'.format(approx_dispersion) + r' $\AA$/pixel') + # ax2.set_ylim(ax1.get_ylim()) + # ax1_ticks = ax1.get_yticks() + # ax2.set_yticklabels([str(round(tick / approx_dispersion, 2)) for tick in ax1_ticks]) + # plt.savefig('{}/interp_vs_our_wls.png'.format(plot_path)) + # plt.close() else: raise ValueError('Only set up to perform Legendre fits currently! Please set fit_type to "Legendre"') From 272c3b84b8a75423deb9fc8ea802046500eda9e3 Mon Sep 17 00:00:00 2001 From: Andrew Howard Date: Mon, 18 Dec 2023 13:26:23 -0800 Subject: [PATCH 006/153] increased NCPU and updated comment --- recipes/quicklook_date.recipe | 4 ++-- scripts/launch_QLP.sh | 2 +- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/recipes/quicklook_date.recipe b/recipes/quicklook_date.recipe index ed3ba9182..4a8437be4 100644 --- a/recipes/quicklook_date.recipe +++ b/recipes/quicklook_date.recipe @@ -1,8 +1,8 @@ # This recipe is run in non-watch mode and is useful for bulk (re)processing # of QLP data products. It computes all QLP data products (from L0, 2D, L1, and L2) # for a given datecode (e.g., 20230711). Data products are put into standard -# locations in /data/QLP///. Bulk reprocessing can be -# done using 'parallel', e.g. a command like > cat