diff --git a/.github/workflows/python-package.yml b/.github/workflows/python-package.yml index 9c8d6090..66e4844d 100644 --- a/.github/workflows/python-package.yml +++ b/.github/workflows/python-package.yml @@ -15,7 +15,7 @@ jobs: runs-on: ubuntu-latest strategy: matrix: - python-version: ['3.8', '3.9', '3.10', '3.11'] + python-version: ['3.10', '3.11', '3.12'] steps: - uses: actions/checkout@v4 @@ -48,5 +48,3 @@ jobs: with: token: ${{ secrets.CODECOV_TOKEN }} # not required for public repos fail_ci_if_error: true # optional (default = false) - - diff --git a/docs/change_history.rst b/docs/change_history.rst index 9f773f93..e91b4909 100644 --- a/docs/change_history.rst +++ b/docs/change_history.rst @@ -1,6 +1,7 @@ Change History ############## + .. _v1.3.8: V1.3.8 Not Released @@ -111,7 +112,7 @@ V1.3.0 06-03-2020 - Removed bias overscan and trimming correction on master bias creation. - Bugs Fixed: - + `--max-targets` was not being used, missed connection in `MainApp`. + + `--max-targets` was not being used, missed connection in `RedSpec`. - Updated keyword values of all reference lamps in the library according to [#292] - Refactored `wavelength.WavelengthCalibration` class moving several methods to diff --git a/docs/conf.py b/docs/conf.py index 0bbc8496..c5e497d4 100644 --- a/docs/conf.py +++ b/docs/conf.py @@ -186,8 +186,6 @@ 'Miscellaneous'), ] - - # -- Options for Epub output ---------------------------------------------- # Bibliographic Dublin Core info. @@ -207,5 +205,3 @@ # A list of files that should not be packed into the epub file. epub_exclude_files = ['search.html'] - - diff --git a/environment.yml b/environment.yml index 091a4347..4dc0a12f 100644 --- a/environment.yml +++ b/environment.yml @@ -14,5 +14,6 @@ dependencies: - pydata-sphinx-theme - astroplan - mock + - setuptools - sphinx - sphinxcontrib.napoleon diff --git a/goodman_pipeline/__init__.py b/goodman_pipeline/__init__.py index b5efb05f..dfd0af4b 100644 --- a/goodman_pipeline/__init__.py +++ b/goodman_pipeline/__init__.py @@ -1,3 +1,5 @@ +# flake8: noqa + from __future__ import absolute_import from .version import __version__ @@ -6,4 +8,4 @@ from . import images from . import core -from .core import setup_logging \ No newline at end of file +from .core import setup_logging diff --git a/goodman_pipeline/core/__init__.py b/goodman_pipeline/core/__init__.py index 32698dc8..eb09d6dc 100644 --- a/goodman_pipeline/core/__init__.py +++ b/goodman_pipeline/core/__init__.py @@ -1,3 +1,5 @@ +# flake8: noqa + from __future__ import (absolute_import, division, print_function, unicode_literals) @@ -13,7 +15,7 @@ # import of functions in core.py from .core import (astroscrappy_lacosmic, - add_wcs_keys, + add_linear_wcs_keys, add_linear_wavelength_solution, bin_reference_data, call_cosmic_rejection, @@ -52,6 +54,7 @@ recenter_broad_lines, recenter_lines, record_trace_information, + record_wavelength_solution_evaluation, save_extracted, search_comp_group, setup_logging, diff --git a/goodman_pipeline/core/core.py b/goodman_pipeline/core/core.py index b9926bd2..1af146a5 100644 --- a/goodman_pipeline/core/core.py +++ b/goodman_pipeline/core/core.py @@ -29,11 +29,13 @@ from astropy.stats import sigma_clip from astropy.time import Time from astroscrappy import detect_cosmics +from ccdproc import CCDData from matplotlib import pyplot as plt from scipy import signal, interpolate from threading import Timer from . import check_version +from ..wcs import WCS __version__ = __import__('goodman_pipeline').__version__ @@ -61,7 +63,7 @@ def astroscrappy_lacosmic(ccd, red_path=None, save_mask=False): return ccd -def add_wcs_keys(ccd): +def add_linear_wcs_keys(ccd): """Adds generic keyword for linear wavelength solution to the header Linear wavelength solutions require a set of standard fits keywords. Later @@ -146,7 +148,7 @@ def add_linear_wavelength_solution(ccd, x_axis, reference_lamp, crpix=1): astronomical tool. (e.g. IRAF) Args: - ccd (CCDData) Instance of :class:`~astropy.nddata.CCDData` + ccd (CCDData): Instance of :class:`~astropy.nddata.CCDData` x_axis (ndarray): Linearized x-axis in angstrom reference_lamp (str): Name of lamp used to get wavelength solution. crpix (int): reference pixel for defining wavelength solution. @@ -162,6 +164,22 @@ def add_linear_wavelength_solution(ccd, x_axis, reference_lamp, crpix=1): new_crpix = crpix new_crval = x_axis[new_crpix - crpix] new_cdelt = x_axis[new_crpix] - x_axis[new_crpix - crpix] + incoming_header_keys = ccd.header.keys() + linear_wcs_keys = ['BANDID1', + 'WCSDIM', + 'CTYPE1', + 'CRVAL1', + 'CRPIX1', + 'CDELT1', + 'CD1_1', + 'LTM1_1', + 'WAT0_001', + 'WAT1_001', + 'DC-FLAG', + 'DCLOG1'] + + if not all([key in incoming_header_keys for key in linear_wcs_keys]): + ccd = add_linear_wcs_keys(ccd=ccd) ccd.header.set('BANDID1', 'spectrum - background none, weights none, ' 'clean no') @@ -175,7 +193,7 @@ def add_linear_wavelength_solution(ccd, x_axis, reference_lamp, crpix=1): ccd.header.set('WAT0_001', 'system=equispec') ccd.header.set('WAT1_001', 'wtype=linear label=Wavelength units=angstroms') ccd.header.set('DC-FLAG', 0) - ccd.header.set('DCLOG1', 'REFSPEC1 = {:s}'.format(reference_lamp)) + ccd.header.set('DCLOG1', 'REFSPEC1 = {:s}'.format(os.path.basename(reference_lamp))) return ccd @@ -323,7 +341,7 @@ def call_cosmic_rejection(ccd, _create(instrument=_instrument, binning=_binning, path=red_path) - #out_prefix = prefix + out_prefix #Move line here + # out_prefix = prefix + out_prefix #Move line here full_path = os.path.join(red_path, f"{out_prefix}_{image_name}") @@ -574,6 +592,7 @@ def cross_correlation(reference, compared, slit_size, serial_binning, + selection_bias='none', mode='full', plot=False): """Do cross correlation of two 1D spectra @@ -591,6 +610,8 @@ def cross_correlation(reference, compared (array): Array to be matched. A new reference lamp. slit_size (float): Slit width in arcseconds serial_binning (int): Binning in the spectral axis + selection_bias (str): For arrays expected to be similar therefore a relatively small shift is expected select 'center'. + 'none' (default) or 'center'. mode (str): Correlation mode for `signal.correlate`. plot (bool): Switch debugging plots on or off. @@ -619,19 +640,35 @@ def cross_correlation(reference, ccorr = signal.correlate(cyaxis1, cyaxis2, mode=mode) - max_index = np.argmax(ccorr) - x_ccorr = np.linspace(-int(len(ccorr) / 2.), int(len(ccorr) / 2.), len(ccorr)) + if selection_bias == 'center': + gaussian_model = models.Gaussian1D(amplitude=1, mean=len(ccorr) / 2., stddev=50) + gaussian_weights = gaussian_model(range(len(x_ccorr))) + gaussian_weighted = gaussian_weights * ccorr + + max_index = np.argmax(gaussian_weighted) + elif selection_bias == 'none': + max_index = np.argmax(ccorr) + + else: + raise NotImplementedError(f"'selection_bias' {selection_bias} is not valid. Options are 'none' and 'center'") correlation_value = x_ccorr[max_index] + + log.debug(f"Found correlation value of {correlation_value}") + if plot: # pragma: no cover plt.ion() plt.title('Cross Correlation') plt.xlabel('Lag Value') plt.ylabel('Correlation Value') - plt.plot(x_ccorr, ccorr) + plt.plot(x_ccorr, ccorr, label='Original Cross Correlation') + if selection_bias == 'center': + plt.plot(x_ccorr, gaussian_weights, label="Centered Gaussian") + plt.plot(x_ccorr, gaussian_weighted, label='Gaussian Weighted') + plt.legend(loc='best') plt.draw() plt.pause(2) plt.clf() @@ -742,9 +779,9 @@ def classify_spectroscopic_data(path, search_pattern): for i in pifc.index.tolist(): radeg, decdeg = ra_dec_to_deg(pifc.obsra.iloc[i], pifc.obsdec.iloc[i]) - pifc.iloc[i, pifc.columns.get_loc('radeg')] = '{:.2f}'.format(radeg) + pifc.iloc[i, pifc.columns.get_loc('radeg')] = '{:.6f}'.format(radeg) - pifc.iloc[i, pifc.columns.get_loc('decdeg')] = '{:.2f}'.format(decdeg) + pifc.iloc[i, pifc.columns.get_loc('decdeg')] = '{:.6f}'.format(decdeg) # now we can compare using degrees confs = pifc.groupby(['slit', @@ -973,7 +1010,7 @@ def kill_process(process): # pragma: no cover # wait for dcr to terminate # dcr.wait() - #go back to the original directory. Could be the same. + # go back to the original directory. Could be the same. os.chdir(cwd) # If no error stderr is an empty string @@ -1016,6 +1053,7 @@ def kill_process(process): # pragma: no cover os.unlink(full_path_out) return ccd + def define_trim_section(sample_image, technique): """Get the initial trim section @@ -1047,9 +1085,7 @@ def define_trim_section(sample_image, technique): # serial binning - dispersion binning # parallel binning - spatial binning spatial_length, dispersion_length = ccd.data.shape - serial_binning, \ - parallel_binning = [int(x) for x - in ccd.header['CCDSUM'].split()] + serial_binning, parallel_binning = [int(x) for x in ccd.header['CCDSUM'].split()] # Trim section is valid for Blue and Red Camera Binning 1x1 and # Spectroscopic ROI @@ -1273,20 +1309,20 @@ def extract_fractional_pixel(ccd, target_trace, target_fwhm, extraction_width, background = background_1 if background_info_1 is None: background_info_1 = "{:.2f}:{:.2f} column {:d}".format( - low_1, high_1, i+1) + low_1, high_1, i + 1) elif background_1 is None and background_2 is not None: background = background_2 if background_info_2 is None: background_info_2 = "{:.2f}:{:.2f} column {:d}".format( - low_2, high_2, i+1) + low_2, high_2, i + 1) else: background = np.mean([background_1, background_2]) if background_info_1 is None: background_info_1 = "{:.2f}:{:.2f} column {:d}".format( - low_1, high_1, i+1) + low_1, high_1, i + 1) if background_info_2 is None: background_info_2 = "{:.2f}:{:.2f} column {:d}".format( - low_2, high_2, i+1) + low_2, high_2, i + 1) # actual background subtraction background_subtracted_column_sum = column_sum - background @@ -1552,8 +1588,11 @@ def get_lines_in_lamp(ccd, plots=False): fig, ax = plt.subplots() fig.canvas.manager.set_window_title('Lines Detected') - mng = plt.get_current_fig_manager() - mng.window.showMaximized() + try: + mng = plt.get_current_fig_manager() + mng.window.showMaximized() + except AttributeError as error: + log.debug(error) ax.set_title('Lines detected in Lamp\n' '{:s}'.format(lamp_header['OBJECT'])) @@ -2071,7 +2110,7 @@ def interpolate_spectrum(spectrum, interpolation_size): """Creates an interpolated version of the input spectrum This method creates an interpolated version of the input array, it is - used mainly for a spectrum but it can also be used with any + used mainly for a spectrum, but it can also be used with any unidimensional array. The reason for doing interpolation is that it allows to find the lines and its respective center more precisely. @@ -2103,7 +2142,7 @@ def is_file_saturated(ccd, threshold): """Detects a saturated image It counts the number of pixels above the saturation_threshold level, then finds - which percentage they represents and if it is above the threshold it + which percentage they represent and if it is above the threshold it will return True. The percentage threshold can be set using the command line argument ``--saturation_threshold``. @@ -2483,6 +2522,7 @@ def read_fits(full_path, technique='Unknown'): GSP_EXTR: Extraction window at first column GSP_BKG1: First background extraction zone GSP_BKG2: Second background extraction zone + GSP_LINE: Data has been linearized. GSP_WRMS: Wavelength solution RMS Error. GSP_WPOI: Number of points used to calculate the wavelength solution Error. @@ -2589,6 +2629,11 @@ def read_fits(full_path, technique='Unknown'): value='none', comment='Second background extraction zone') + if 'GSP_LINE' not in all_keys: + ccd.header.set('GSP_LINE', + value='FALSE', + comment='Data has been linearized') + if 'GSP_WRMS' not in all_keys: ccd.header.set('GSP_WRMS', value='none', @@ -2817,6 +2862,26 @@ def record_trace_information(ccd, trace_info): return ccd +def record_wavelength_solution_evaluation(ccd: CCDData, rms_error: float, n_points: int, n_rejections: int): + """Add record of wavalength solution evaluation + + Args: + ccd (CCDData): data to be updated. + rms_error (float): Root Nean Square Error of the wavelength solution. + n_points (int): Number of points used to evaluate the solution. + n_rejections (int): Number of points rejected while evaluating the solution. + + Returns: + ccd modified + + """ + + ccd.header.set('GSP_WRMS', value=rms_error) + ccd.header.set('GSP_WPOI', value=n_points) + ccd.header.set('GSP_WREJ', value=n_rejections) + return ccd + + def save_extracted(ccd, destination, prefix='e', target_number=1): """Save extracted spectrum while adding a prefix. @@ -3745,6 +3810,9 @@ def __init__(self, reference_dir): reference_dir (str): full path to the reference data directory """ self.log = logging.getLogger(__name__) + if not os.path.isdir(reference_dir): + self.log.error(f"Given `reference_dir`: {reference_dir} does not exist.") + self.log.debug(f"Current Directory: {os.getcwd()}") self.reference_dir = reference_dir reference_collection = ccdproc.ImageFileCollection(self.reference_dir) self.ref_lamp_collection = reference_collection.summary.to_pandas() @@ -3764,6 +3832,689 @@ def __init__(self, reference_dir): 'LAMP_DOM', 'LAMP_DPE'] + self.lamps_elements = {'LAMP_HGA': ['Hg', 'Ar'], + 'LAMP_NE': ['Ne'], + 'LAMP_AR': ['Ar'], + 'LAMP_FE': ['Fe', 'He'], + 'LAMP_CU': ['Cu', 'He']} + + self.line_list_files = {'cu': 'Cu_3000A-10000A_clean.csv', + 'he': 'He_3000A-10000A_clean.csv', + 'ne': 'Ne_3000A-10000A_clean.csv', + 'ar': 'Ar_3000A-10000A_clean.csv', + 'hg': 'Hg_3000A-10000A_clean.csv'} + + self.line_list = {'hg': [3125.67, + 3131.70, + 3341.48, + 3650.153, + 3654.84, + 3663.279, + 4046.563, + 4077.831, + 4358.327, + 4916.068, + 5460.750, + 5769.598, + 5790.663], + 'cu': [1999.698, + 2035.854, + 2037.127, + 2043.802, + 2135.981, + 2165.09, + 2178.94, + 2179.410, + 2181.72, + 2189.630, + 2192.268, + 2199.58, + 2210.268, + 2214.58, + 2218.108, + 2225.70, + 2227.78, + 2230.08, + 2242.618, + 2244.26, + 2247.002, + 2263.08, + 2293.84, + 2392.63, + 2406.66, + 2441.64, + 2492.15, + 2618.37, + 2689.300, + 2700.962, + 2703.184, + 2713.508, + 2718.778, + 2766.37, + 2769.669, + 2877.100, + 2961.16, + 3036.10, + 3063.41, + 3108.60, + 3194.10, + 3243.16, + 3247.54, + 3273.96, + 3290.54, + 3307.95, + 3337.84, + 3530.38, + 3686.555, + 4062.64, + 4651.12, + 4909.734, + 4931.698, + 4953.724, + 5051.793, + 5105.54, + 5153.24, + 5218.20, + 5292.52, + 5700.24, + 5782.13, + 6000.120, + 6150.384, + 6154.222, + 6216.939, + 6219.844, + 6273.349, + 6301.009, + 6377.840, + 6423.884, + 6448.559, + 6470.168, + 6481.437, + 6624.292, + 6641.396, + 7404.354, + 7652.333, + 7664.648, + 7778.738, + 7805.184, + 7807.659, + 7825.654, + 7902.553, + 7933.13, + 7988.163, + 8092.63, + 8283.160, + 8511.061, + 9861.280, + 9864.137, + 10054.938], + 'ne': [5852.488, + 5881.895, + 5944.834, + 5975.534, + 6029.997, + 6074.337, + 6096.163, + 6128.450, + 6143.062, + 6163.594, + 6217.281, + 6266.495, + 6304.789, + 6334.428, + 6382.991, + 6402.246, + 6506.528, + 6532.882, + 6598.953, + 6652.092, + 6678.276, + 6717.043, + 6929.467, + 7024.050, + 7032.413, + 7173.938, + 7245.166, + 7438.898, + 7488.871, + 7535.774, + 7544.044, + 7943.180, + 8082.457, + 8136.406, + 8300.324, + 8377.606, + 8418.426, + 8495.359, + 8591.258, + 8634.647, + 8654.384, + 8853.867, + 8865.756, + 8919.499, + 8988.58, + 9148.68, + 9201.76, + 9300.85, + 9313.98, + 9326.52, + 9354.218, + 9425.38, + 9459.21, + 9486.680, + 9534.167, + 9547.40, + 9665.424], + 'he': [3187.745, + 3354.55, + 3613.643, + 3819.6072, + 3888.648, + 3964.7289, + 4009.268, + 4120.815, + 4168.967, + 4437.551, + 4471.479, + 4713.1455, + 4921.931, + 5015.6779, + 5047.738, + 5875.621, + 7281.349, + 10830.337], + 'ar': [3093.4019, + 3139.0176, + 3161.3726, + 3169.6685, + 3181.0376, + 3204.3210, + 3243.6887, + 3249.8003, + 3281.7016, + 3307.2283, + 3350.9243, + 3376.4359, + 3388.5309, + 3454.0952, + 3476.7474, + 3478.2324, + 3480.5055, + 3499.4765, + 3509.7785, + 3514.3877, + 3519.9936, + 3521.2601, + 3521.9781, + 3535.3196, + 3548.5144, + 3554.3058, + 3556.9041, + 3559.5081, + 3561.0304, + 3565.0298, + 3567.6564, + 3576.6156, + 3581.6084, + 3582.3546, + 3588.4407, + 3605.8792, + 3606.5218, + 3622.1375, + 3632.6831, + 3637.0310, + 3639.8329, + 3649.8323, + 3650.8896, + 3655.2782, + 3656.0498, + 3659.5289, + 3660.437, + 3669.6024, + 3671.0048, + 3673.2645, + 3678.2701, + 3680.0609, + 3682.5448, + 3690.8951, + 3706.9302, + 3709.9088, + 3714.7337, + 3717.1713, + 3718.2065, + 3720.4265, + 3724.5165, + 3729.3087, + 3737.889, + 3746.4476, + 3750.4799, + 3753.5177, + 3754.0498, + 3763.5053, + 3765.27, + 3766.1186, + 3770.52, + 3780.8398, + 3786.3824, + 3796.5934, + 3799.382, + 3803.1724, + 3808.5748, + 3809.4561, + 3819.0159, + 3825.6729, + 3826.8072, + 3834.6787, + 3841.5187, + 3844.7311, + 3845.4055, + 3850.5813, + 3856.1382, + 3861.747, + 3868.5284, + 3872.1371, + 3875.2645, + 3880.3332, + 3891.4017, + 3891.9792, + 3894.66, + 3895.2502, + 3900.6266, + 3911.576, + 3914.7675, + 3925.7188, + 3926.0456, + 3928.6233, + 3931.2359, + 3932.5466, + 3944.2717, + 3946.0971, + 3947.5046, + 3948.9789, + 3952.7291, + 3958.38, + 3968.3594, + 3974.4766, + 3974.759, + 3979.3559, + 3979.7155, + 3988.1576, + 3992.0535, + 3994.7918, + 3999.252, + 4001.1379, + 4005.3628, + 4011.2093, + 4013.8566, + 4019.8429, + 4022.629, + 4031.3783, + 4033.8093, + 4035.46, + 4038.8043, + 4042.8937, + 4044.4179, + 4045.6773, + 4045.9654, + 4047.4815, + 4052.9208, + 4054.5258, + 4062.641, + 4063.238, + 4070.7835, + 4072.0047, + 4072.3849, + 4076.6284, + 4076.9432, + 4079.5738, + 4082.3872, + 4099.4563, + 4103.9121, + 4112.8153, + 4116.3753, + 4128.64, + 4129.6823, + 4131.7235, + 4144.2435, + 4156.086, + 4158.5905, + 4164.1795, + 4168.9682, + 4178.3658, + 4179.2973, + 4181.8836, + 4189.6511, + 4190.7129, + 4191.0294, + 4198.317, + 4199.8891, + 4200.6745, + 4201.5549, + 4201.9715, + 4203.4109, + 4217.4308, + 4218.6649, + 4222.6373, + 4226.6089, + 4226.9876, + 4228.158, + 4229.8696, + 4237.2198, + 4248.956, + 4251.1846, + 4255.6034, + 4259.3619, + 4266.2864, + 4266.5271, + 4272.1689, + 4277.5282, + 4282.8976, + 4297.9645, + 4300.1008, + 4300.6495, + 4309.2392, + 4331.1995, + 4332.0297, + 4333.5612, + 4335.3379, + 4337.0708, + 4338.2314, + 4345.168, + 4345.8966, + 4348.064, + 4352.2049, + 4362.0662, + 4363.7945, + 4367.8316, + 4370.7532, + 4371.329, + 4372.49, + 4374.8579, + 4375.9542, + 4379.6668, + 4382.928, + 4383.7535, + 4385.0566, + 4386.9656, + 4397.7971, + 4400.0968, + 4400.9863, + 4401.755, + 4404.9022, + 4406.4704, + 4420.9124, + 4423.9944, + 4426.0011, + 4430.189, + 4430.9963, + 4433.838, + 4438.1175, + 4439.4614, + 4439.8793, + 4440.1216, + 4445.8483, + 4448.4597, + 4448.8792, + 4449.5206, + 4460.5574, + 4474.7594, + 4480.35, + 4481.8107, + 4490.9816, + 4498.5384, + 4502.9268, + 4507.8339, + 4509.374, + 4510.7332, + 4522.323, + 4530.5523, + 4530.785, + 4535.4903, + 4537.6426, + 4543.8692, + 4545.0519, + 4547.7589, + 4561.0128, + 4563.7429, + 4564.4054, + 4579.3495, + 4589.8978, + 4596.0967, + 4598.7627, + 4609.5673, + 4628.4409, + 4637.2328, + 4651.124, + 4657.9012, + 4682.2759, + 4702.3161, + 4721.5910, + 4726.8683, + 4732.0532, + 4735.9058, + 4764.8646, + 4806.0205, + 4847.8095, + 4865.9105, + 4876.2611, + 4879.8635, + 4882.2432, + 4888.2612, + 4889.0422, + 4904.7516, + 4914.3146, + 4933.2091, + 4942.9214, + 4965.0795, + 4972.1597, + 5009.3344, + 5017.1628, + 5062.0371, + 5090.4951, + 5118.2023, + 5125.7654, + 5141.7827, + 5145.3083, + 5151.3907, + 5162.2846, + 5165.7728, + 5176.2292, + 5187.7462, + 5216.8139, + 5218.2020, + 5221.2710, + 5254.4648, + 5421.3517, + 5439.9891, + 5442.2427, + 5443.6893, + 5451.6520, + 5457.4157, + 5467.1608, + 5473.4516, + 5490.1194, + 5495.8738, + 5498.1841, + 5506.1128, + 5524.9570, + 5558.7020, + 5572.5413, + 5577.6845, + 5581.8714, + 5588.7200, + 5597.4756, + 5601.1216, + 5606.7330, + 5641.3751, + 5648.6863, + 5650.7043, + 5659.1272, + 5681.9001, + 5691.6612, + 5700.8730, + 5738.3869, + 5739.5196, + 5772.1143, + 5774.0087, + 5783.5360, + 5802.0798, + 5834.2633, + 5860.3103, + 5882.6242, + 5888.5841, + 5912.0853, + 5916.5992, + 5927.1258, + 5928.8130, + 5942.6686, + 5949.2583, + 5971.6008, + 5998.9987, + 6005.7242, + 6013.6777, + 6025.1500, + 6032.1274, + 6043.2233, + 6046.8977, + 6052.7229, + 6059.3725, + 6064.7508, + 6081.2433, + 6085.8797, + 6090.7848, + 6098.8031, + 6101.1615, + 6103.5390, + 6105.6351, + 6113.4657, + 6114.9234, + 6119.6565, + 6123.3619, + 6127.4160, + 6128.7227, + 6145.4411, + 6155.2385, + 6165.1232, + 6170.1740, + 6173.0964, + 6174.4032, + 6212.5031, + 6215.9383, + 6243.1201, + 6296.8722, + 6307.6570, + 6333.1459, + 6357.0229, + 6364.8937, + 6369.5748, + 6384.7169, + 6396.6097, + 6399.2065, + 6403.0128, + 6416.3071, + 6418.3703, + 6431.5550, + 6437.6003, + 6441.8994, + 6443.8598, + 6466.5526, + 6468.0483, + 6472.4294, + 6481.1453, + 6483.0825, + 6493.9694, + 6499.1061, + 6538.1120, + 6604.8534, + 6620.9665, + 6632.0837, + 6638.2207, + 6639.7403, + 6643.6976, + 6656.9386, + 6660.6761, + 6664.0510, + 6666.3588, + 6677.2817, + 6684.2929, + 6719.2184, + 6752.8335, + 6756.1631, + 6766.6117, + 6827.2488, + 6861.2688, + 6863.5350, + 6871.2891, + 6879.5824, + 6888.1742, + 6937.6642, + 6951.4776, + 6965.4307, + 7030.2514, + 7067.2181, + 7147.0416, + 7206.9804, + 7272.9359, + 7311.7159, + 7353.2930, + 7372.1184, + 7383.9805, + 7503.8691, + 7514.6518, + 7635.1060, + 7723.7611, + 7724.2072, + 7948.1764, + 8006.1567, + 8014.7857, + 8103.6931, + 8115.311, + 8264.5225, + 8408.2096, + 8424.6475, + 8521.4422, + 8667.9442, + 9122.9674, + 9224.4992, + 9354.2198, + 9657.7863, + 9784.5028, + 10470.0535]} + + def get_line_list_by_lamp(self, ccd): + """Get the reference lines for elements in the lamp's name + + Splits the name in chunks of two characters assuming each one of them + represents an element in the comparison lamp, then fetches the list of + line positions available and appends it to a list that will be return + ordered. + + Args: + lamp_name (str): Lamp's name as in the header keyword OBJECT. + + Returns: + line_list(list): Sorted line list. + + """ + + + if ccd.header['OBSTYPE'] in ['ARC', 'COMP']: + elements = [] + for lamp_key in ccd.header['LAMP_*']: + print(lamp_key, ccd.header[lamp_key]) + if ccd.header[lamp_key] == 'TRUE': + elements.extend(self.lamps_elements[lamp_key]) + # elements = [lamp_name[i:i + 2].lower() for i in range(0, + # len(lamp_name), + # 2)] + print(elements) + line_list = [] + for element in elements: + line_list.extend(self.line_list[element.lower()]) + return sorted(line_list) + else: + log.critical(f"Obstype {ccd.header['OBSTYPE']} it's not expected to have lamp signal.") + def get_reference_lamp(self, header): """Finds a suitable template lamp from the catalog @@ -3772,7 +4523,7 @@ def get_reference_lamp(self, header): lamp. Returns: - full path to best matching reference lamp. + full path to best matching reference lamp or None. """ @@ -3786,7 +4537,7 @@ def get_reference_lamp(self, header): (self.ref_lamp_collection['lamp_cu'] == header['LAMP_CU']) & (self.ref_lamp_collection['wavmode'] == header['WAVMODE']))] if filtered_collection.empty: - error_message = "Unable to find a match for: "\ + error_message = "Unable to find a match in the reference library for: "\ "LAMP_HGA = {}, "\ "LAMP_NE = {}, "\ "LAMP_AR = {}, "\ @@ -3799,7 +4550,7 @@ def get_reference_lamp(self, header): header['LAMP_CU'], header['WAVMODE']) self.log.error(error_message) - raise NoMatchFound(error_message) + return None else: filtered_collection = self.ref_lamp_collection[ (self.ref_lamp_collection['object'] == header['object']) & @@ -3812,7 +4563,7 @@ def get_reference_lamp(self, header): "WAVMODE = {}".format(header['OBJECT'], header['WAVMODE']) self.log.error(error_message) - raise NoMatchFound(error_message) + return None if len(filtered_collection) == 1: self.log.info( @@ -3919,6 +4670,54 @@ def check_comp_group(self, comp_group): return comp_group return None + def create_reference_lamp(self, ccd, wavelength_solution, lines_pixel, lines_angstrom): + self.log.info(f"Creating Reference lamp.") + assert len(lines_pixel) == len(lines_angstrom), f"Pixel and Angstrom list must have same length." + self.lines_pixel = lines_pixel + self.lines_angstrom = lines_angstrom + + wcs = WCS() + + ccd = self._record_lines(ccd=ccd) + ccd = wcs.write_gsp_wcs(ccd=ccd, model=wavelength_solution) + + lamp_elements = [] + for keyword in ccd.header['LAMP_*']: + if ccd.header[keyword] == 'TRUE': + lamp_elements.extend(self.lamps_elements[keyword]) + + elements = set(lamp_elements) + + lamp_elements_name = "".join(elements) + + if ccd.header['FILTER2'] in ['NO_FILTER']: + pass + + new_name = f"goodman_comp_{ccd.header['WAVMODE']}_{ccd.header['FILTER2']}_{lamp_elements_name}.fits" + + self.log.info(f"Created new reference lamp with name {new_name}") + + write_fits(ccd=ccd, full_path=new_name, overwrite=True) + + def _record_lines(self, ccd): + self.log.info("Adding pixel and angstrom line record.") + last_added = 'GSP_WREJ' + for i in range(len(self.lines_pixel)): + keyword_name = f"GSP_P{i + 1:03}" + ccd.header.set(keyword_name, + value=self.lines_pixel[i], + comment="Line location in pixel value", + after=last_added) + last_added = keyword_name + for i in range(len(self.lines_angstrom)): + keyword_name = f"GSP_A{i + 1:03}" + ccd.header.set(keyword_name, + value=self.lines_angstrom[i], + comment="Line location in angstrom value", + after=last_added) + last_added = keyword_name + return ccd + def _recover_lines(self): """Read lines from the reference lamp's header.""" self.log.info("Recovering line information from reference Lamp.") @@ -3931,6 +4730,7 @@ def _recover_lines(self): assert pixel_key[-3:] == angstrom_key[-3:] assert angstrom_key in self._ccd.header if int(float(self._ccd.header[angstrom_key])) != 0: + self.log.debug(f"Retrieving Line {re.sub('GSP_P', '', pixel_key)}") self.lines_pixel.append(float(self._ccd.header[pixel_key])) self.lines_angstrom.append( float(self._ccd.header[angstrom_key])) @@ -4715,10 +5515,10 @@ def _fit_gaussian(fitter, profile_model = [] if stddev_min is None: stddev_min = 0 - log.debug(f"Setting STDDEV minimum value to {stddev_min} pixels. Set it with `--profile-min-width`.") + log.debug(f"Setting STDDEV minimum value to {stddev_min} pixels. Set it with `--target-min-width`.") if stddev_max is None: stddev_max = 4 * order - log.debug(f"Setting STDDEV maximum value to {stddev_max} pixels. Set it with `--profile-max-width`.") + log.debug(f"Setting STDDEV maximum value to {stddev_max} pixels. Set it with `--target-max-width`.") log.debug(f"Using minimum STDDEV = {stddev_min} pixels.") log.debug(f"Using maximum STDDEV = {stddev_max} pixels.") for peak in selected_peaks: diff --git a/goodman_pipeline/core/tests/test_core.py b/goodman_pipeline/core/tests/test_core.py index fcefdc60..f611c4a7 100644 --- a/goodman_pipeline/core/tests/test_core.py +++ b/goodman_pipeline/core/tests/test_core.py @@ -34,7 +34,7 @@ # import of functions in core.py from ..core import (astroscrappy_lacosmic, add_linear_wavelength_solution, - add_wcs_keys, + add_linear_wcs_keys, bias_subtract, bin_reference_data, call_cosmic_rejection, @@ -93,7 +93,7 @@ def setUp(self): self.ccd = CCDData(data=np.random.random_sample(200), meta=fits.Header(), unit='adu') - self.ccd = add_wcs_keys(ccd=self.ccd) + self.ccd = add_linear_wcs_keys(ccd=self.ccd) self.ccd.header.set('SLIT', value='1.0_LONG_SLIT', comment="slit [arcsec]") @@ -145,9 +145,7 @@ def test_add_wcs_keys(self): 'DC-FLAG', 'DCLOG1'] - - - self.test_ccd = add_wcs_keys(ccd=self.test_ccd) + self.test_ccd = add_linear_wcs_keys(ccd=self.test_ccd) for key in wcs_keys: self.assertIn(key, self.test_ccd.header) @@ -208,8 +206,6 @@ def test__bin_reference_data(self): self.assertEqual(len(new_wavelength), np.floor(len(wavelength) / i)) - - class CentralWavelength(TestCase): def setUp(self): @@ -1840,9 +1836,9 @@ def test_get_reference_lamp_no_match_status_keys(self): self.ccd.header.set('WAVMODE', value='400_M2') - self.assertRaises(NoMatchFound, - self.rd.get_reference_lamp, - self.ccd.header) + reference_lamp = self.rd.get_reference_lamp(header=self.ccd.header) + + self.assertIsNone(reference_lamp) def test_get_reference_lamp_exist_with_object_key(self): self.ccd.header.set('OBJECT', value='HgArNe') @@ -1868,9 +1864,9 @@ def test_get_reference_lamp_does_not_exist(self): self.ccd.header.set('OBJECT', value='HgArCu') self.ccd.header.set('WAVMODE', value='400_M5') - self.assertRaises(NoMatchFound, - self.rd.get_reference_lamp, - self.ccd.header) + reference_lamp = self.rd.get_reference_lamp(header=self.ccd.header) + + self.assertIsNone(reference_lamp) def test_lamp_exist(self): self.ccd.header.set('LAMP_HGA', value='TRUE') diff --git a/goodman_pipeline/images/__init__.py b/goodman_pipeline/images/__init__.py index 20bff635..286a3bbc 100644 --- a/goodman_pipeline/images/__init__.py +++ b/goodman_pipeline/images/__init__.py @@ -1,6 +1,4 @@ -"""Goodman CCD Reduction Tool +# flake8: noqa -""" - -from .goodman_ccd import MainApp, get_args +from .goodman_ccd import RedCCD, get_args from .image_processor import ImageProcessor diff --git a/goodman_pipeline/images/data_classifier.py b/goodman_pipeline/images/data_classifier.py index acf22eac..089c92fc 100644 --- a/goodman_pipeline/images/data_classifier.py +++ b/goodman_pipeline/images/data_classifier.py @@ -143,18 +143,18 @@ def _get_obs_technique(self): wavmodes = [str(w).upper() for w in self.objects_collection.wavmode.unique()] if len(wavmodes) == 1 and wavmodes[0] == 'IMAGING': - self.technique = 'Imaging' + self.technique = 'Imaging' elif 'IMAGING' in wavmodes and len(wavmodes) > 1: - self.log.error('There seems to be Imaging and Spectroscopic ' - 'data. I will assume the Imaging data are ' - 'acquisition images therefore they will be ' - 'ignored.') - self.log.info("If you really have Imaging data, please process " - "them in a separated folder.") - self.technique = 'Spectroscopy' + self.log.error('There seems to be Imaging and Spectroscopic ' + 'data. I will assume the Imaging data are ' + 'acquisition images therefore they will be ' + 'ignored.') + self.log.info("If you really have Imaging data, please process " + "them in a separated folder.") + self.technique = 'Spectroscopy' else: - self.technique = 'Spectroscopy' + self.technique = 'Spectroscopy' # inform the results, no need to return self.log.info('Detected {:s} Data from {:s} ' 'Camera'.format(self.technique, self.instrument)) diff --git a/goodman_pipeline/images/goodman_ccd.py b/goodman_pipeline/images/goodman_ccd.py index 3072ddca..dc2b9670 100755 --- a/goodman_pipeline/images/goodman_ccd.py +++ b/goodman_pipeline/images/goodman_ccd.py @@ -152,10 +152,10 @@ def get_args(arguments=None): return args -class MainApp(object): +class RedCCD(object): def __init__(self): - """This method initializes the MainApp class + """This method initializes the RedCCD class The main task of this method is to call the get_args function that returns an argparse object. @@ -174,7 +174,7 @@ def __init__(self): self._pipeline_version = __version__ def __call__(self, args=None): - """Call method for MainApp + """Call method for RedCCD From the arguments this method finds the raw_path attribute and checks its contents for the existence of files containing the '.fits' string. @@ -392,5 +392,5 @@ def _check_args(self): if __name__ == '__main__': # pragma: no cover - main_app = MainApp() - main_app() + RED_CCD = RedCCD() + RED_CCD() diff --git a/goodman_pipeline/images/night_organizer.py b/goodman_pipeline/images/night_organizer.py index 3f34af0f..8354a51b 100644 --- a/goodman_pipeline/images/night_organizer.py +++ b/goodman_pipeline/images/night_organizer.py @@ -95,7 +95,7 @@ def __call__(self): _imaging_file = self.file_collection[ self.file_collection.wavmode == 'IMAGING'] if not _imaging_file.empty: - self.log.warning("Ignoring all Imaging data. Assuming they are " + self.log.warning("Ignoring all files where `WAVMODE` is `IMAGING`. Assuming they are " "not science exposures.") self.log.info("If any of the files listed below is a science " "file, please process it/them in a separate " @@ -108,8 +108,8 @@ def __call__(self): self.file_collection.wavmode != 'IMAGING'].reset_index(drop=True) elif self.technique == 'Imaging': - self.log.warning("Ignoring all files where `wavmode` is not " - "Imaging.") + self.log.warning("Ignoring all files where `WAVMODE` is not " + "IMAGING.") self.file_collection = self.file_collection[ self.file_collection.wavmode == 'IMAGING'].reset_index(drop=True) diff --git a/goodman_pipeline/images/tests/test_goodman_ccd.py b/goodman_pipeline/images/tests/test_goodman_ccd.py index 91d62e1f..e677e5c0 100644 --- a/goodman_pipeline/images/tests/test_goodman_ccd.py +++ b/goodman_pipeline/images/tests/test_goodman_ccd.py @@ -1,18 +1,18 @@ from __future__ import absolute_import from unittest import TestCase, skip -from ..goodman_ccd import get_args, MainApp +from ..goodman_ccd import get_args, RedCCD class MainAppTest(TestCase): def setUp(self): - self.main_app = MainApp() + self.red_ccd = RedCCD() def test___call__(self): - self.assertRaises(SystemExit, self.main_app) + self.assertRaises(SystemExit, self.red_ccd) def test___call___show_version(self): arguments = ['--version'] args = get_args(arguments=arguments) - self.assertRaises(SystemExit, self.main_app, args) + self.assertRaises(SystemExit, self.red_ccd, args) diff --git a/goodman_pipeline/scripts/redccd b/goodman_pipeline/scripts/redccd index 4f1115ca..b2ce6d9f 100755 --- a/goodman_pipeline/scripts/redccd +++ b/goodman_pipeline/scripts/redccd @@ -3,7 +3,7 @@ from __future__ import (absolute_import, division, print_function, unicode_literals) from goodman_pipeline.core import setup_logging -from goodman_pipeline.images import MainApp +from goodman_pipeline.images import RedCCD import sys @@ -13,5 +13,5 @@ if '-h' not in sys.argv and \ setup_logging() if __name__ == '__main__': - GOODMAN_CCD = MainApp() + GOODMAN_CCD = RedCCD() GOODMAN_CCD() diff --git a/goodman_pipeline/scripts/redspec b/goodman_pipeline/scripts/redspec index dfed4911..028cd2a5 100755 --- a/goodman_pipeline/scripts/redspec +++ b/goodman_pipeline/scripts/redspec @@ -3,7 +3,7 @@ from __future__ import (absolute_import, division, print_function, unicode_literals) from goodman_pipeline.core import setup_logging -from goodman_pipeline.spectroscopy import MainApp +from goodman_pipeline.spectroscopy import RedSpec import sys @@ -14,5 +14,5 @@ if '-h' not in sys.argv and \ if __name__ == '__main__': - GOODMAN_SPECTROSCOPY = MainApp() + GOODMAN_SPECTROSCOPY = RedSpec() GOODMAN_SPECTROSCOPY() diff --git a/goodman_pipeline/spectroscopy/__init__.py b/goodman_pipeline/spectroscopy/__init__.py index 3674e444..1ce52103 100644 --- a/goodman_pipeline/spectroscopy/__init__.py +++ b/goodman_pipeline/spectroscopy/__init__.py @@ -1,6 +1,6 @@ -"""Goodman Spectroscopic Tools +# flake8: noqa -""" from __future__ import absolute_import -from .redspec import MainApp, get_args +from .redspec import RedSpec, get_args from .wavelength import WavelengthCalibration +from ._interactive_wavelength import InteractiveWavelengthCalibration diff --git a/goodman_pipeline/spectroscopy/_interactive_wavelength.py b/goodman_pipeline/spectroscopy/_interactive_wavelength.py new file mode 100644 index 00000000..2c25bc7f --- /dev/null +++ b/goodman_pipeline/spectroscopy/_interactive_wavelength.py @@ -0,0 +1,1577 @@ +from __future__ import (absolute_import, division, print_function, + unicode_literals) + +import logging +import os +import re +import sys + +import astropy.units as u +import matplotlib.pyplot as plt +import matplotlib.ticker as mtick +import numpy as np +from astropy.convolution import convolve, Gaussian1DKernel +from astropy.modeling import models, fitting +from astropy.stats import sigma_clip +from ccdproc import CCDData + +from ..wcs.wcs import WCS +from ..core import (evaluate_wavelength_solution, + get_lines_in_lamp, + get_spectral_characteristics, + record_wavelength_solution_evaluation) + +from ..core import ReferenceData + +FORMAT = '%(levelname)s:%(filename)s:%(module)s: %(message)s' +logging.basicConfig(level=logging.DEBUG, format=FORMAT) +log = logging.getLogger(__name__) + + +SHOW_PLOTS = False + + +class InteractiveWavelengthCalibration(object): + """Wavelength Calibration Class + + The WavelengthCalibration class is instantiated for each of the science + images, which are treated as a "science object". In this first release it + can find a wavelength solution for a given comparison lamp using an + interactive GUI based on Matplotlib. Although it works very good, for the + next release there is a plan for creating an independent GUI based on QT in + order to work better in different screen sizes and other topic such as + showing warnings, messages and help. + + This class takes 1D spectrum with no wavelength calibration and returns fits + files with wavelength solutions using the FITS standard for linear + solutions. Goodman spectra are slightly non-linear therefore they are + linearized and smoothed before they are returned for the user. + + """ + + def __init__(self): + """Wavelength Calibration Class Initialization + + A WavelengthCalibration class is instantiated for each science target + being processed, i.e. every science image. + + Notes: + This class violates some conventions as for length and number of + attributes is concerned. Solving this is part of a prioritary plans + for next release. + + Args: + args (object): Runtime arguments. + + """ + + # TODO - Documentation missing + self.lamp = None + self.poly_order = 2 + self.wcs = WCS() + self._wsolution = None + self._rms_error = None + self._n_rejections = None + self._n_points = None + # print(self.args.reference_dir) + self.reference_data = None + # self.science_object = science_object + self.slit_offset = None + self.interpolation_size = 200 + self.line_search_method = 'derivative' + """Instrument configuration and spectral characteristics""" + self.pixel_size = 15 * u.micrometer + self.pixel_scale = 0.15 * u.arcsec + self.goodman_focal_length = 377.3 * u.mm + self.grating_frequency = None + self.grating_angle = None + self.camera_angle = None + self.serial_binning = None + self.parallel_binning = None + + self.pixel_count = None + self.alpha = None + self.beta = None + self.center_wavelength = None + self.blue_limit = None + self.red_limit = None + """Interactive wavelength finding""" + self.reference_marks_x = [] + self.reference_marks_y = [] + self.raw_data_marks_x = [] + self.raw_data_marks_y = [] + self.click_input_enabled = True + self.reference_bb = None + self.raw_data_bb = None + self.contextual_bb = None + self.i_fig = None + self.ax1 = None + self.ax2 = None + self.ax3 = None + self.ax4 = None + self.ax4_plots = None + self.ax4_com = None + self.ax4_rlv = None + self.legends = None + self.points_ref = None + self.points_raw = None + self.line_raw = None + self.line_list = [] + self.line_pixels = [] + self.line_angstroms = [] + self.ref_filling_value = 1000 + self.raw_filling_value = 1000 + self.events = True + self.first = True + self.evaluation_comment = None + # self.binning = self.lamp.header[] + self.pixel_center = [] + """this data must come parsed""" + # self.science_pack = sci_pack + # self.sci_filename = self.science_object.file_name + # self.history_of_lamps_solutions = {} + self._reference_lamp = None + self.reference_solution = None + + def __call__(self, + ccd: CCDData, + comparison_lamp: CCDData, + save_data_to, + reference_data, + object_number=None, + output_prefix='w', + wsolution_obj=None, + linearize=False, + plot_results=False, + save_plots=False, + plots=False): + """Call method for the WavelengthSolution Class + + It takes extracted data and produces wavelength calibrated 1D FITS file. + The call method takes care of the order and logic needed to call the + different methods. A wavelength solution can be recycled for the next + science object. In that case, the wavelength solution is parsed as an + argument and then there is no need to calculate it again. The recycling + part has to be implemented in the caller function. + + Args: + ccd (object): a ccdproc.CCDData instance + comp_list (list): Comparison lamps for the science target that will + be processed here. Every element of this list is an instance of + ccdproc.CCDData. + object_number (int): In case of multiple detections in a single + image this number will be added as a suffix before `.fits` in + order to allow for multiple 1D files. Default value is None. + wsolution_obj (object): Mathematical model of the wavelength + solution if exists. If it doesn't is a None + + Returns: + wavelength_solution (object): The mathematical model of the + wavelength solution. If it fails to create it will return a + None element. + + """ + assert isinstance(ccd, CCDData) + assert isinstance(comparison_lamp, CCDData) + + self._save_data_to = save_data_to + + if os.path.isdir(reference_data): + self.reference_data = ReferenceData(reference_dir=reference_data) + else: + log.error(f"Reference Data Directory provided :{reference_data} does not exist.") + + self.i_fig = None + + log.info('Processing Science Target: ' + '{:s}'.format(ccd.header['OBJECT'])) + + self.lamp = comparison_lamp + + self.line_list = self.reference_data.get_line_list_by_lamp(ccd=self.lamp) + try: + self.calibration_lamp = self.lamp.header['GSP_FNAM'] + except KeyError: + self.calibration_lamp = '' + + self.raw_pixel_axis = range(self.lamp.shape[0]) + self.lamp_name = self.lamp.header['OBJECT'] + + log.info('Processing Comparison Lamp: ' + '{:s}'.format(self.lamp_name)) + self.lines_center = get_lines_in_lamp(ccd=self.lamp) + self.spectral = get_spectral_characteristics( + ccd=self.lamp, + pixel_size=self.pixel_size, + instrument_focal_length=self.goodman_focal_length) + object_name = ccd.header['OBJECT'] + # try: + self._interactive_wavelength_solution( + object_name=object_name) + # except TypeError as error: + # log.error(error) + + if self._wsolution is not None: + if plots or plot_results or save_plots: + + plt.close(1) + if not plots: + plt.ion() + # plt.show() + else: + plt.ioff() + + wavelength_axis = self._wsolution(range(ccd.data.size)) + + object_name = ccd.header['OBJECT'] + grating = ccd.header['GRATING'] + + fig_title = 'Wavelength Calibrated Data : ' \ + '{:s}\n{:s}'.format(object_name, grating) + + fig, ax1 = plt.subplots(1) + fig.canvas.set_window_title(ccd.header['GSP_FNAM']) + # ax1 = fig.add_subplot(111) + manager = plt.get_current_fig_manager() + if plt.get_backend() == u'GTK3Agg': + manager.window.maximize() + elif plt.get_backend() == u'Qt5Agg': + manager.window.showMaximized() + + ax1.set_title(fig_title) + ax1.set_xlabel('Wavelength (Angstrom)') + ax1.set_ylabel('Intensity (ADU)') + ax1.set_xlim((wavelength_axis[0], wavelength_axis[-1])) + # plt.close(1) + + ax1.plot(wavelength_axis, + ccd.data, + color='k', + label='Data') + + ax1.legend(loc='best') + fig.tight_layout() + if save_plots: + log.info('Saving plots') + plots_dir = os.path.join(save_data_to, 'plots') + if not os.path.isdir(plots_dir): + os.mkdir(plots_dir) + plot_name = re.sub('.fits', + '.png', + ccd.header['GSP_FNAM']) + plot_path = os.path.join(plots_dir, plot_name) + # print(plot_path) + plt.savefig(plot_path, dpi=300) + log.info('Saved plot as {:s} file ' + 'DPI=300'.format(plot_name)) + + if plot_results: + plt.show() + else: + plt.draw() + plt.pause(2) + plt.ioff() + + return self._wsolution + else: + log.error('It was not possible to get a wavelength ' + 'solution from this lamp.') + return None + + @property + def reference_lamp(self): + return self._reference_lamp + + @property + def wavelength_solution(self): + return self._wsolution + + @property + def wavelength_solution_evaluation(self): + return self._rms_error, self._n_points, self._n_rejections + + def create_reference_lamp(self): + log.info("Saving as template") + self.lamp = record_wavelength_solution_evaluation(ccd=self.lamp, + rms_error=self._rms_error, + n_points=self._n_points, + n_rejections=self._n_rejections) + + self.reference_data.create_reference_lamp( + ccd=self.lamp, + wavelength_solution=self._wsolution, + lines_pixel=self.line_pixels, + lines_angstrom=self.line_angstroms) + + def get_best_filling_value(self, data): + """Find the best y-value to locate marks + + The autmatically added points will be placed at a fixed location in the + y-axis. This value is calculated by doing a 2-sigma clipping with 5 + iterations. Then the masked out values are removed and the median is + calculated. + + Args: + data (array): Array of 1D data + + Returns: + Median value of clipped data. + + """ + clipped_data = sigma_clip(data, sigma=2, maxiters=5) + clean_data = clipped_data[~clipped_data.mask] + log.debug("Found best filling value" + " at {:f}".format(np.median(clean_data))) + return np.median(clean_data) + + def recenter_lines(self, data, lines, plots=False): + """Finds the centroid of an emission line + + For every line center (pixel value) it will scan left first until the + data stops decreasing, it assumes it is an emission line and then will + scan right until it stops decreasing too. Defined those limits it will + use the line data in between and calculate the centroid. + + Notes: + This method is used to recenter relatively narrow lines only, there + is a special method for dealing with broad lines. + + Args: + data (array): numpy.ndarray instance. or the data attribute of a + ccdproc.CCDData instance. + lines (list): A line list in pixel values. + plots (bool): If True will plot spectral line as well as the input + center and the recentered value. + + Returns: + A list containing the recentered line positions. + + """ + new_center = [] + x_size = data.shape[0] + median = np.median(data) + for line in lines: + # TODO (simon): Check if this definition is valid, so far is not + # TODO (cont..): critical + left_limit = 0 + right_limit = 1 + condition = True + left_index = int(line) + + while condition and left_index - 2 > 0: + + if (data[left_index - 1] > data[left_index]) and \ + (data[left_index - 2] > data[left_index - 1]): + + condition = False + left_limit = left_index + + elif data[left_index] < median: + condition = False + left_limit = left_index + + else: + left_limit = left_index + + left_index -= 1 + + # id right limit + condition = True + right_index = int(line) + while condition and right_index + 2 < x_size - 1: + + if (data[right_index + 1] > data[right_index]) and \ + (data[right_index + 2] > data[right_index + 1]): + + condition = False + right_limit = right_index + + elif data[right_index] < median: + condition = False + right_limit = right_index + + else: + right_limit = right_index + right_index += 1 + index_diff = [abs(line - left_index), abs(line - right_index)] + + sub_x_axis = range(line - min(index_diff), + (line + min(index_diff)) + 1) + + sub_data = data[line - min(index_diff):(line + min(index_diff)) + 1] + centroid = np.sum(sub_x_axis * sub_data) / np.sum(sub_data) + + # checks for asymmetries + differences = [abs(data[line] - data[left_limit]), + abs(data[line] - data[right_limit])] + + if max(differences) / min(differences) >= 2.: + if plots: + plt.axvspan(line - 1, line + 1, color='g', alpha=0.3) + new_center.append(line) + else: + new_center.append(centroid) + if plots: + fig = plt.figure(9) + fig.canvas.set_window_title('Lines Detected in Lamp') + plt.axhline(median, color='b') + + plt.plot(self.raw_pixel_axis, + data, + color='k', + label='Lamp Data') + + for line in lines: + + plt.axvline(line + 1, + color='k', + linestyle=':', + label='First Detected Center') + + for center in new_center: + + plt.axvline(center, + color='k', + linestyle='.-', + label='New Center') + + plt.show() + return new_center + + @staticmethod + def recenter_broad_lines(lamp_data, lines, order): + """Recenter broad lines + + Notes: + This method is used to recenter broad lines only, there is a special + method for dealing with narrower lines. + + Args: + lamp_data (array): numpy.ndarray instance. It contains the lamp + data. + lines (list): A line list in pixel values. + order (float): A rough estimate of the FWHM of the lines in pixels + in the data. It is calculated using the slit size divided by the + pixel scale multiplied by the binning. + + Returns: + A list containing the recentered line positions. + + """ + # TODO (simon): use slit size information for a square function + # TODO (simon): convolution + new_line_centers = [] + gaussian_kernel = Gaussian1DKernel(stddev=2.) + lamp_data = convolve(lamp_data, gaussian_kernel) + for line in lines: + lower_index = max(0, int(line - order)) + upper_index = min(len(lamp_data), int(line + order)) + lamp_sample = lamp_data[lower_index:upper_index] + x_axis = np.linspace(lower_index, upper_index, len(lamp_sample)) + line_max = np.max(lamp_sample) + + gaussian_model = models.Gaussian1D(amplitude=line_max, + mean=line, + stddev=order) + + fit_gaussian = fitting.LevMarLSQFitter() + fitted_gaussian = fit_gaussian(gaussian_model, x_axis, lamp_sample) + new_line_centers.append(fitted_gaussian.mean.value) + # if self.args.debug_mode: + # plt.plot(x_axis, lamp_sample) + # plt.plot(x_axis, gaussian_model(x_axis)) + # plt.plot(x_axis, fitted_gaussian(x_axis), color='k') + # plt.axvline(line) + # plt.show() + return new_line_centers + + def recenter_line_by_data(self, data_name, x_data): + """Finds a better center for a click-selected line + + This method is called by another method that handles click events. An + argument is parsed that will tell which plot was clicked and what is + the x-value in data coordinates. Then the closest pixel center will be + found and from there will extract a 20 pixel wide sample of the data + (this could be a future improvement: the width of the extraction should + depend on the FWHM of the lines). The sample of the data is used to + calculate a centroid (center of mass) which is a good approximation but + could be influenced by data shape or if the click was too far + (unquantified yet). That is why for the reference data, a database of + laboratory line center will be + used and for raw data the line centers are calculated earlier in the + process, independently of any human input. + + It wil also plot the sample, the centroid and the reference line center + at the fourth subplot, bottom right corner. + + Args: + data_name (str): 'reference' or 'raw-data' is where the click was + done + x_data (float): click x-axis value in data coordinates + + Returns: + reference_line_value (float): The value of the line center that will + be used later to do the wavelength fit + + """ + if data_name == 'reference': + print(f"{self.reference_solution=}") + + reference_line_index = np.argmin( + abs(self.line_list - x_data)) + + reference_line_value = self.line_list[reference_line_index] + + if not self.reference_solution: + return reference_line_value + + pseudo_center = np.argmin(abs(self.reference_solution[0] - x_data)) + + sub_x = self.reference_solution[0][ + pseudo_center - 10: pseudo_center + 10] + + sub_y = self.reference_solution[1][ + pseudo_center - 10: pseudo_center + 10] + + center_of_mass = np.sum(sub_x * sub_y) / np.sum(sub_y) + # print 'centroid ', center_of_mass + # plt.figure(3) + # if self.ax4_plots is not None or self.ax4_com is not None or + # self.ax4_rlv is not None: + try: + self.ax4.cla() + self.ax4.relim() + except NameError as err: + log.error(err) + + self.ax4.set_title('Reference Data Clicked Line') + self.ax4.set_xlabel('Wavelength (Angstrom)') + self.ax4.set_ylabel('Intensity (Counts)') + + self.ax4_plots = self.ax4.plot(sub_x, + sub_y, + color='k', + label='Data') + + self.ax4_rlv = self.ax4.axvline(reference_line_value, + linestyle='-', + color='r', + label='Reference Line Value') + + self.ax4_com = self.ax4.axvline(center_of_mass, + linestyle='--', + color='b', + label='Centroid') + + self.ax4.yaxis.set_major_formatter(mtick.FormatStrFormatter('%.1e')) + self.ax4.legend(loc=3, framealpha=0.5) + self.i_fig.canvas.draw() + # return center_of_mass + return reference_line_value + elif data_name == 'raw-data': + pseudo_center = np.argmin(abs(self.raw_pixel_axis - x_data)) + raw_line_index = np.argmin(abs(self.lines_center - x_data)) + raw_line_value = self.lines_center[raw_line_index] + # print(raw_line_value) + sub_x = self.raw_pixel_axis[pseudo_center - 10: pseudo_center + 10] + sub_y = self.lamp.data[pseudo_center - 10: pseudo_center + 10] + center_of_mass = np.sum(sub_x * sub_y) / np.sum(sub_y) + # print 'centroid ', center_of_mass + # plt.figure(3) + # if self.ax4_plots is not None or self.ax4_com is not None or + # self.ax4_rlv is not None: + try: + self.ax4.cla() + self.ax4.relim() + except NameError as err: + log.error(err) + self.ax4.set_title('Raw Data Clicked Line') + self.ax4.set_xlabel('Pixel Axis') + self.ax4.set_ylabel('Intensity (Counts)') + + self.ax4_plots = self.ax4.plot(sub_x, + sub_y, + color='k', + label='Data') + + self.ax4_rlv = self.ax4.axvline(raw_line_value, + linestyle='-', + color='r', + label='Line Center') + self.ax4_com = self.ax4.axvline(center_of_mass, + linestyle='--', + color='b', + label='Centroid') + + self.ax4.yaxis.set_major_formatter(mtick.FormatStrFormatter('%.1e')) + self.ax4.legend(loc=3, framealpha=0.5) + self.i_fig.canvas.draw() + # return center_of_mass + return raw_line_value + else: + log.error('Unrecognized data name') + + def _interactive_wavelength_solution(self, object_name=''): + """Find the wavelength solution interactively + + This method uses the graphical capabilities of matplotlib in particular + the user interface (UI) capabilities such as click and key-pressed + events. We could say that it implements a matplotlib Graphical User + Interface. There are some limitation though, for instance, the access to + a better formatted help text. + + It will display four plots of the same height but their width have a + 4:1 ratio being the leftmost the wider. + + Top Left Plot: Displays the spectrum of the lamp we want to calibrate, + or as it is called in the plot, the Raw Data, therefore the axis are + intensity (ADU) versus dispersion (pixels). There will be red vertical + dashed lines that are lines detected by pipeline. If you see the a line + without a red dashed line it means it was not detected. + + Bottom Left Plot: Displays the reference lamp data, therefore the axis + are intensity (ADU) versus wavelength (Angstrom). In this case there are + vertical dotted lines plotted with a relatively low alpha value. They + represent reference laboratory lines obtained from the NIST database at + https://physics.nist.gov/PhysRefData/ASD/lines_form.html. + + Top Right Plot: This is a quick reference help text. The formatting of + the text is quite difficult here due to matplotlib's own limitations. + This text is also displayed on the terminal. + + Bottom Right Plot: This is a dynamic plot since it shows different kind + of information. First it shows short messages, such as warnings, errors, + or general information. Once you select a line it will show a zoomed + version of the line plus the reference value chosen as well as were the + click was done. And finally, it shows the dispersion of the wavelength + solution. This could be one of the most useful plots in the entire + screen. + + Usage: A quick reference of the instructions is shown in the Top Right + Plot, but in general you have to select the matching lines in both data + plots this will create a table of pixels versus angstrom which at its + roots a wavelength solution. Then this will have to be translated to a + mathematical model. Having this mathematical model the wavelength + solution can be applied, evaluated etc. Among the features implemented + you can: + - Delete a line or a matching pair. + - Find more lines automatically once you have a decent wavelength + solution. + - Remove all the automatically added lines at once, this is useful in + some cases when adding more lines actually makes the solution worst. + - Find the RMS Error of the solution + - Remove ALL the points at a given point to start over. + - Fit a wavelength solution + + + + Notes: + This method uses the Qt5Agg backend, in theory it could also work + with GTK3Agg but it is being forced to use Qt5Agg. + + """ + log.debug("Starting interactive wavelength calibration") + plt.switch_backend('Qt5Agg') + + # disable full screen to allow the use of f for fitting the solution + + plt.rcParams['keymap.fullscreen'] = [u'ctrl+f'] + + try: + self._reference_lamp = self.reference_data.get_reference_lamp( + header=self.lamp.header) + except NotImplementedError: + self._reference_lamp = self.reference_data.get_reference_lamps_by_name( + lamp_name=self.lamp.header['OBJECT']) + log.warning('Could not find a perfect match for reference ' + 'data') + # reference_file = None + # log.critical('Could not find a comparison lamp in the + # reference.') + + # reference_file = self.reference_data.get_reference_lamps_by_name( + # self.lamp_name) + + if self._reference_lamp is not None: + log.info('Using reference file: {:s}'.format(self._reference_lamp.header['GSP_FNAM'])) + reference_plots_enabled = True + + self.reference_solution = self.wcs.read_gsp_wcs(ccd=self._reference_lamp) + print(self.reference_solution) + else: + reference_plots_enabled = False + log.error('Please Check the OBJECT Keyword of your reference ' + 'data') + + # update filling value + self.raw_filling_value = self.get_best_filling_value( + data=self.lamp.data) + + if self._reference_lamp is not None: + self.ref_filling_value = self.get_best_filling_value(data=self._reference_lamp.data) + + # ------- Plots ------- + self.i_fig, ((self.ax1, self.ax2), (self.ax3, self.ax4)) = \ + plt.subplots(2, + 2, + gridspec_kw={'width_ratios': [4, 1]}) + + self.i_fig.canvas.setWindowTitle('Science Target: {:s}'.format( + object_name)) + + manager = plt.get_current_fig_manager() + if plt.get_backend() == u'GTK3Agg': + manager.window.maximize() + elif plt.get_backend() == u'Qt5Agg': + manager.window.showMaximized() + + self.ax1.set_title('Raw Data - {:s}\n{:s} - {:s}'.format( + self.lamp_name, + self.lamp.header['GRATING'], + self.lamp.header['SLIT'])) + + self.ax1.set_xlabel('Pixels') + self.ax1.set_ylabel('Intensity (counts)') + self.ax1.plot([], linestyle='--', color='r', label='Detected Lines') + for idline in self.lines_center: + self.ax1.axvline(idline, linestyle='--', color='r') + + self.ax1.plot(self.raw_pixel_axis, + self.lamp.data, + color='k', + label='Raw Data') + + self.ax1.set_xlim((0, len(self.lamp.data))) + self.ax1.yaxis.set_major_formatter(mtick.FormatStrFormatter('%.1e')) + self.ax1.legend(loc=2) + + # Update y limits to have an extra 5% to top and bottom + ax1_ylim = self.ax1.get_ylim() + ax1_y_range = ax1_ylim[1] - ax1_ylim[0] + + self.ax1.set_ylim((ax1_ylim[0] - 0.05 * ax1_y_range, + ax1_ylim[1] + 0.05 * ax1_y_range)) + + self.ax3.set_title('Reference Data - {:s}'.format(self.lamp_name)) + self.ax3.set_xlabel('Wavelength (Angstrom)') + self.ax3.set_ylabel('Intensity (counts)') + self.ax3.set_xlim((self.spectral['blue'].value, self.spectral['red'].value)) + self.ax3.yaxis.set_major_formatter(mtick.FormatStrFormatter('%.1e')) + + self.ax3.plot([], + linestyle=':', + color='r', + label='Reference Line Values') + + for rline in self.line_list: + self.ax3.axvline(rline, linestyle=':', color='r') + + if reference_plots_enabled: + + self.ax3.plot(self.reference_solution[0], + self.reference_solution[1], + color='k', + label='Reference Lamp Data') + + self.ax3.legend(loc=2) + + # Update y limits to have an extra 5% to top and bottom + ax3_ylim = self.ax3.get_ylim() + ax3_y_range = ax3_ylim[1] - ax3_ylim[0] + + self.ax3.set_ylim((ax3_ylim[0] - 0.05 * ax3_y_range, + ax3_ylim[1] + 0.05 * ax3_y_range)) + # print(ax3_ylim) + + self.display_help_text() + + # zoomed plot + self.display_onscreen_message('Use middle button click to select a ' + 'line') + + plt.subplots_adjust(left=0.05, + right=0.99, + top=0.96, + bottom=0.04, + hspace=0.17, + wspace=0.11) + + self.raw_data_bb = self.ax1.get_position() + self.reference_bb = self.ax3.get_position() + self.contextual_bb = self.ax4.get_position() + + # if self.click_input_enabled: + self.i_fig.canvas.mpl_connect('button_press_event', self._on_click) + self.i_fig.canvas.mpl_connect('key_press_event', self._key_pressed) + # print self._wsolution + plt.show() + return True + + def _on_click(self, event): + """Handles Click events for Interactive Mode + + Calls the method register_mark + + Args: + event (object): Click event + """ + if event.button == 3: + self._register_mark(event) + # TODO (simon): Make sure the text below is useless + # else: + # print(event.button) + # elif event.button == 3: + # if len(self.reference_marks_x) == len(self.raw_data_marks_x): + # self.click_input_enabled = False + # log.info('Leaving interactive mode') + # else: + # if len(self.reference_marks_x) < len(self.raw_data_marks_x): + # log.info('There is {:d} click missing in the ' + # 'Reference plot'.format( + # len(self.raw_data_marks_x) - + # len(self.reference_marks_x))) + # else: + # log.info('There is {:d} click missing in the ' + # 'New Data plot'.format( + # len(self.reference_marks_x) - + # len(self.raw_data_marks_x))) + + def _key_pressed(self, event): + """Key event handler + + There are several key events that need to be taken care of. + See a brief description of each one of them below: + + F1 or ?: Prints a help message + F2 or f: Fit wavelength solution model. + F3 or a: Find new lines. + F4: Evaluate solution + F6 or l: Linearize data although this is done automatically after the + wavelength function is fit + d: deletes closest point + ctrl+d: deletes all recorded marks + ctrl+z: Reverts the action of F3 or a. + Middle Button Click or m: records data location. + Enter: Close figure and apply solution if exists. + Shift+Enter: Close the program with sys.exit(0) + + Notes: + This method must be simplified + + Args: + event (object): Key pressed event + + """ + self.events = True + if event.key == 'f1' or event.key == '?': + log.info('Print help regarding interactive mode') + print("F1 or ?: Prints Help.") + print("F2 or f: Fit wavelength solution model.") + print("F3 or a: Find new lines.") + print("F4: Evaluate solution") + print("F6 or l: Linearize data (for testing not definitive)") + print("d: deletes closest point") + # print("l : resample spectrum to a linear dispersion axis") + print("ctrl+d: deletes all recorded marks") + print("ctrl+t: Save as template") + print("ctrl+z: Go back to previous solution " + "(deletes automatic added points") + print('Right Button Click or p: records data location.') + print("Enter: Close figure and apply solution if exists.") + elif event.key == 'f2' or event.key == 'f': + log.debug('Calling function to fit wavelength Solution') + self.fit_pixel_to_wavelength() + self._plot_raw_over_reference() + elif event.key == 'f3' or event.key == 'a': + if self._wsolution is not None: + self._find_more_lines() + self._update_marks_plot('reference') + self._update_marks_plot('raw_data') + else: + log.debug('Wavelength solution is None') + elif event.key == 'f4': + if self._wsolution is not None and len(self.raw_data_marks_x) > 0: + self.evaluate_solution(plots=True) + elif event.key == 'f5' or event.key == 'd': + # TODO (simon): simplify this code. + + figure_x, figure_y = \ + self.i_fig.transFigure.inverted().transform((event.x, event.y)) + + if self.raw_data_bb.contains(figure_x, figure_y): + log.debug('Deleting raw point') + # print abs(self.raw_data_marks_x - event.xdata) a[:] = + + closer_index = int(np.argmin( + [abs(list_val - event.xdata) for list_val in + self.raw_data_marks_x])) + + # print 'Index ', closer_index + if len(self.raw_data_marks_x) == len(self.reference_marks_x): + self.raw_data_marks_x.pop(closer_index) + self.raw_data_marks_y.pop(closer_index) + self.reference_marks_x.pop(closer_index) + self.reference_marks_y.pop(closer_index) + self._update_marks_plot('reference') + self._update_marks_plot('raw_data') + else: + if closer_index == len(self.raw_data_marks_x) - 1: + self.raw_data_marks_x.pop(closer_index) + self.raw_data_marks_y.pop(closer_index) + self._update_marks_plot('raw_data') + else: + self.raw_data_marks_x.pop(closer_index) + self.raw_data_marks_y.pop(closer_index) + self.reference_marks_x.pop(closer_index) + self.reference_marks_y.pop(closer_index) + self._update_marks_plot('reference') + self._update_marks_plot('raw_data') + + elif self.reference_bb.contains(figure_x, figure_y): + log.debug('Deleting reference point') + # print 'reference ', self.reference_marks_x, self.re + # print self.reference_marks_x + # print abs(self.reference_marks_x - event.xdata) + + closer_index = int(np.argmin( + [abs(list_val - event.xdata) for list_val in + self.reference_marks_x])) + + if len(self.raw_data_marks_x) == len(self.reference_marks_x): + self.raw_data_marks_x.pop(closer_index) + self.raw_data_marks_y.pop(closer_index) + self.reference_marks_x.pop(closer_index) + self.reference_marks_y.pop(closer_index) + self._update_marks_plot('reference') + self._update_marks_plot('raw_data') + else: + if closer_index == len(self.reference_marks_x) - 1: + self.reference_marks_x.pop(closer_index) + self.reference_marks_y.pop(closer_index) + self._update_marks_plot('reference') + else: + self.raw_data_marks_x.pop(closer_index) + self.raw_data_marks_y.pop(closer_index) + self.reference_marks_x.pop(closer_index) + self.reference_marks_y.pop(closer_index) + self._update_marks_plot('reference') + self._update_marks_plot('raw_data') + + elif self.contextual_bb.contains(figure_x, figure_y): + log.warning("Can't delete points from here because points " + "represent each detected line in the raw data.") + # closer_index_ref = int(np.argmin( + # [abs(list_val - event.ydata) for list_val in + # self.reference_marks_x])) + # + # closer_index_raw = int(np.argmin( + # [abs(list_val - event.xdata) for list_val in + # self.raw_data_marks_x])) + # + # log.debug('Raw Index {:d}, Ref Index {:d}'.format( + # closer_index_raw, + # closer_index_ref)) + # self.raw_data_marks_x.pop(closer_index_raw) + # self.raw_data_marks_y.pop(closer_index_raw) + # self.reference_marks_x.pop(closer_index_raw) + # self.reference_marks_y.pop(closer_index_raw) + # self.update_marks_plot('reference') + # self.update_marks_plot('raw_data') + # self.update_marks_plot('eval_plots') + # elif event.key == 'i': + # figure_x, figure_y = self.i_fig.transFigure.inverted().transform( + # (event.x, event.y)) + # if self.contextual_bb.contains(figure_x, figure_y): + # log.debug("Trying to identify point.") + # + # closer_x_index = int(np.argmin( + # [abs(list_val - event.xdata) for list_val in + # self.raw_data_marks_x])) + # print(self.raw_data_marks_x[closer_x_index]) + # print(self.reference_marks_x[closer_x_index]) + + elif event.key == 'f6' or event.key == 'l': + log.info('Linearize and smoothing spectrum') + if self._wsolution is not None: + self.linearize_spectrum(self.lamp.data, plots=True) + + elif event.key == 'ctrl+z': + log.info('Deleting automatic added points. If exist.') + + if self.raw_data_marks_x is not [] and \ + self.reference_marks_x is not []: + + to_remove = [] + for i in range(len(self.raw_data_marks_x)): + # print self.raw_data_marks[i], self.filling_value + if self.raw_data_marks_y[i] == self.raw_filling_value: + to_remove.append(i) + # print to_remove + to_remove = np.array(sorted(to_remove, reverse=True)) + if len(to_remove) > 0: + for index in to_remove: + self.raw_data_marks_x.pop(index) + self.raw_data_marks_y.pop(index) + self.reference_marks_x.pop(index) + self.reference_marks_y.pop(index) + self._update_marks_plot('reference') + self._update_marks_plot('raw_data') + # else: + # print self.raw_click_plot, self.ref_click_plot, 'mmm' + + elif event.key == 'ctrl+d': + try: + log.info('Deleting all recording Clicks') + self.display_onscreen_message( + message='All points deleted') + self.reference_marks_x = [] + self.reference_marks_y = [] + self.raw_data_marks_x = [] + self.raw_data_marks_y = [] + self._update_marks_plot('delete') + self._plot_raw_over_reference(remove=True) + log.info('All points deleted!') + except: + log.error('No points deleted') + elif event.key == 'ctrl+t': + self.create_reference_lamp() + + elif event.key == 'p': + self._register_mark(event) + + elif event.key == 'enter': + if self._wsolution is not None: + log.info('Closing figure') + plt.close('all') + else: + message = 'There is still no wavelength solution!' + log.info(message) + self.display_onscreen_message(message) + + elif event.key == 'm': + self._register_mark(event) + + elif event.key in ['ctrl+q', 'ctrl+c']: + log.info('Pressed Ctrl+q. Closing the program') + sys.exit(0) + + else: + log.debug("No action for key pressed: {:s}".format(event.key)) + pass + + def _register_mark(self, event): + """Marks a line + + Detects where the click was done or m-key was pressed and calls the + corresponding method. It handles the middle button click and m-key being + pressed. There are two regions of interest as for where a click was + done. The raw and reference data respectively. For any of such regions + it will call the method that recenter the line and once the desired + value is returned it will be appended to the list that contains all the + correspondent line positions, raw (pixels) and reference (angstrom) + + Args: + event (object): Click or m-key pressed event + """ + if event.xdata is not None and event.ydata is not None: + figure_x, figure_y = \ + self.i_fig.transFigure.inverted().transform((event.x, event.y)) + + if self.reference_bb.contains(figure_x, figure_y): + # self.reference_marks.append([event.xdata, event.ydata]) + self.reference_marks_x.append( + self.recenter_line_by_data('reference', event.xdata)) + + self.reference_marks_y.append(event.ydata) + self._update_marks_plot('reference') + elif self.raw_data_bb.contains(figure_x, figure_y): + # self.raw_data_marks.append([event.xdata, event.ydata]) + self.raw_data_marks_x.append( + self.recenter_line_by_data('raw-data', event.xdata)) + + self.raw_data_marks_y.append(event.ydata) + self._update_marks_plot('raw_data') + else: + log.debug('{:f} {:f} Are not contained'.format(figure_x, + figure_y)) + else: + log.error('Clicked Region is out of boundaries') + + def _find_more_lines(self): + """Method to add more lines given that a wavelength solution already + exists + + This method is part of the interactive wavelength solution mechanism. + If a wavelength solution exist it uses the line centers in pixels to + estimate their respective wavelength and then search for the closest + value in the list of reference lines for the elements in the comparison + lamp. Then it filters the worst of them by doing sigma clipping. + Finally it adds them to the class' attributes that contains the list of + reference points. + + Better results are obtained if the solution is already decent. Visual + inspection also improves final result. + """ + new_physical = [] + new_wavelength = [] + square_differences = [] + if self._wsolution is not None: + wlines = self._wsolution(self.lines_center) + for i in range(len(wlines)): + # [abs(list_val - wlines[i]) for list_val in \ + # self.reference_data.get_line_list_by_name(self.lamp_name)] + + closer_index = np.argmin( + abs(self.line_list - wlines[i])) + + rline = self.line_list[closer_index] + + rw_difference = wlines[i] - rline + # print('Difference w - r ', rw_difference, rline) + square_differences.append(rw_difference ** 2) + new_physical.append(self.lines_center[i]) + new_wavelength.append(rline) + clipped_differences = sigma_clip(square_differences, + sigma=2, + maxiters=3) + + if len(new_wavelength) == len(new_physical) == \ + len(clipped_differences): + + for i in range(len(new_wavelength)): + if clipped_differences[i] is not \ + np.ma.masked and new_wavelength[i] not in \ + self.reference_marks_x: + + self.reference_marks_x.append(new_wavelength[i]) + self.reference_marks_y.append(self.ref_filling_value) + self.raw_data_marks_x.append(new_physical[i]) + self.raw_data_marks_y.append(self.raw_filling_value) + return True + + def _update_marks_plot(self, action=None): + """Update the points that represent marks on lamp plots + + When you mark a line a red dot marks the position of the line at the + exact y location of the click, for the x location it will use the value + obtained by means of the recentering method. There are three possible + actions: Update the reference plot's click, the raw data marks or + delete them all. + + Args: + action (str): A string that could be 'reference', 'raw_data' or + 'delete' depending on the action desired + """ + # print(type(action), type(pixel_axis), type(differences)) + if action == 'reference': + if self.points_ref is not None: + try: + self.points_ref.remove() + self.ax3.relim() + log.debug('Removing reference marks') + except: + log.debug('Reference points is None') + pass + log.debug("Plot new marks") + self.points_ref, = self.ax3.plot(self.reference_marks_x, + self.reference_marks_y, + linestyle='None', + marker='o', + color='r') + self.i_fig.canvas.draw() + elif action == 'raw_data': + # print self.points_raw + # print dir(self.points_raw) + if self.points_raw is not None: + try: + self.points_raw.remove() + self.ax1.relim() + except ValueError as err: + log.error(err) + self.points_raw, = self.ax1.plot(self.raw_data_marks_x, + self.raw_data_marks_y, + linestyle='None', + marker='o', + color='r') + self.i_fig.canvas.draw() + elif action == 'delete': + if self.points_raw is not None and self.points_ref is not None: + self.points_raw.remove() + self.ax1.relim() + self.points_ref.remove() + self.ax3.relim() + self.i_fig.canvas.draw() + else: + log.error('Unknown Action {:s}'.format(action)) + + def _plot_raw_over_reference(self, remove=False): + """Overplot raw data over reference lamp using current wavelength + solution model + + Once the wavelength solution is obtained this method is called to apply + the already mentioned solution to the raw data and then overplot it on + the reference lamp plot. This is very useful to have a visual idea of + how far or close the solution is. + + Args: + remove (bool): True or False depending whether you want to remove + the overplotted lamp or not + """ + if self._wsolution is not None: + if self.line_raw is not None: + try: + self.line_raw.remove() + self.ax3.relim() + except: + pass + if not remove: + # TODO(simon): catch TypeError Exception and correct what is + # TODO (cont): causing it + self.line_raw, = self.ax3.plot( + self._wsolution(self.raw_pixel_axis), + self.lamp.data, + linestyle='-', + color='r', + alpha=0.4, + label='New Wavelength Solution') + + self.ax3.legend(loc=2) + self.i_fig.canvas.draw() + + def evaluate_solution(self, plots=False): + """Calculate the Root Mean Square Error of the solution + + Once the wavelength solution is obtained it has to be evaluated. The + line centers found for the raw comparison lamp will be converted to, + according to the new solution, angstrom. Then for each line the closest + reference line value is obtained. The difference is stored. Then these + differences are cleaned by means of a sigma clipping method that will + rule out any outlier or any line that is not well-matched. Then, using + the sigma clipped differences the Root Mean Square error is calculated. + + It also creates a plot in the bottom right subplot of the interactive + window, showing a scatter plot plus some information regarding the + quality of the fit. + + Args: + plots (bool): Whether to create the plot or not + + Returns: + results (list): Contains three elements: _rms_error (float), + npoints (int), _n_rejections (int) + + """ + if self._wsolution is not None: + differences = np.array([]) + wavelength_line_centers = self._wsolution(self.lines_center) + + for wline in wavelength_line_centers: + closer_index = np.argmin( + abs(self.line_list - wline)) + + rline = self.line_list[closer_index] + + rw_difference = wline - rline + differences = np.append(differences, rw_difference) + + clipping_sigma = 2. + + clipped_differences = differences + once_clipped_differences = differences + + old_rms_error = None + if self._rms_error is not None: + old_rms_error = float(self._rms_error) + + self._rms_error, self._n_points, self._n_rejections = evaluate_wavelength_solution( + clipped_differences=clipped_differences) + + log.info('RMS Error : {:.3f}'.format(self._rms_error)) + + if plots: + if self.ax4_plots is not None or \ + self.ax4_com is not None or \ + self.ax4_rlv is not None: + + try: + self.ax4.cla() + self.ax4.relim() + except NameError as err: + log.error(err) + + self.ax4.set_title('RMS Error {:.3f} \n' + '{:d} points ({:d} ' + 'rejected)'.format(self._rms_error, + self._n_points, + self._n_rejections)) + + self.ax4.set_ylim(once_clipped_differences.min(), + once_clipped_differences.max()) + + self.ax4.set_xlim(np.min(self.lines_center), + np.max(self.lines_center)) + + self.ax4_rlv = self.ax4.scatter(self.lines_center, + differences, + marker='x', + color='k', + label='Rejected Points') + + self.ax4_com = self.ax4.axhspan(clipped_differences.min(), + clipped_differences.max(), + color='k', + alpha=0.4, + edgecolor=None, + label='{:.1f} Sigma'.format( + clipping_sigma)) + + self.ax4_plots = self.ax4.scatter(self.lines_center, + clipped_differences, + label='Differences') + + if self._rms_error is not None and old_rms_error is not None: + rms_error_difference = self._rms_error - old_rms_error + + if rms_error_difference > 0.001: + increment_color = 'red' + elif rms_error_difference < -0.001: + increment_color = 'green' + else: + increment_color = 'white' + + message = r'$\Delta$ RMSE {:+.3f}'.format( + rms_error_difference) + + self.ax4.text(0.05, 0.95, + message, + verticalalignment='top', + horizontalalignment='left', + transform=self.ax4.transAxes, + color=increment_color, + fontsize=15) + + self.ax4.set_xlabel('Pixel Axis (Pixels)') + self.ax4.set_ylabel('Difference (Angstroms)') + + self.ax4.legend(loc=3, framealpha=0.5) + self.i_fig.canvas.draw() + + results = [self._rms_error, self._n_points, self._n_rejections] + return results + else: + log.error('Solution is still non-existent!') + + def fit_pixel_to_wavelength(self): + """Does the fit to find the wavelength solution + + Once you have four data points on each side (raw and reference or pixel + and angstrom) it calculates the fit using a Chebyshev model of third + degree. This was chosen because it worked better compared to the rest. + There is a slight deviation from linearity in all Goodman data, + therefore a linear model could not be used, also is said that a Spline + of third degree is "too flexible" which I also experienced and since the + deviation from linearity is not extreme it seemed that it was not + necessary to implement. + + This method checks that the data that will be used as input to calculate + the fit have the same dimensions and warns the user in case is not. + + Returns: + None (None): An empty return is created to finish the execution of + the method when a fit will not be possible + + """ + if len(self.reference_marks_x) and len(self.raw_data_marks_x) > 0: + + if len(self.reference_marks_x) < 4 or \ + len(self.raw_data_marks_x) < 4: + + message = 'Not enough marks! Minimum 4 each side.' + self.display_onscreen_message(message) + return + + if len(self.reference_marks_x) != len(self.raw_data_marks_x): + if len(self.reference_marks_x) < len(self.raw_data_marks_x): + n = len(self.raw_data_marks_x) - len(self.reference_marks_x) + if n == 1: + message_text = '{:d} Reference Click is ' \ + 'missing!.'.format(n) + else: + message_text = '{:d} Reference Clicks are ' \ + 'missing!.'.format(n) + else: + n = len(self.reference_marks_x) - len(self.raw_data_marks_x) + if n == 1: + message_text = '{:d} Raw Click is missing!.'.format(n) + else: + message_text = '{:d} Raw Clicks are missing!.'.format(n) + self.display_onscreen_message(message_text) + else: + self.line_pixels = [] + self.line_angstroms = [] + for i in range(len(self.reference_marks_x)): + self.line_pixels.append(self.raw_data_marks_x[i]) + self.line_angstroms.append(self.reference_marks_x[i]) + + self._wsolution = self.wcs.fit(physical=self.line_pixels, + wavelength=self.line_angstroms, + model_name='chebyshev', + degree=self.poly_order) + + self.evaluate_solution(plots=True) + + else: + log.error('Clicks record is empty') + self.display_onscreen_message(message='Clicks record is empty') + if self._wsolution is not None: + self._wsolution = None + + def add_gsp_wcs(self): + pass + + def display_onscreen_message(self, message='', color='red'): + """Uses the fourth subplot to display a message + + Displays a warning message on the bottom right subplot of the + interactive window. It is capable to break down the message in more + than one line if necessary. + + Args: + message (str): The message to be displayed + color (str): Color name for the font's color + + """ + full_message = [message] + if len(message) > 30: + full_message = [] + split_message = message.split(' ') + line_length = 0 + # new_line = '' + e = 0 + for i in range(len(split_message)): + # print(i, len(split_message)) + line_length += len(split_message[i]) + 1 + if line_length >= 30: + new_line = ' '.join(split_message[e:i]) + # print(new_line, len(new_line)) + full_message.append(new_line) + # new_line = '' + line_length = 0 + e = i + if i == len(split_message) - 1: + new_line = ' '.join(split_message[e:]) + # print(new_line, len(new_line)) + full_message.append(new_line) + + self.ax4.cla() + self.ax4.relim() + self.ax4.set_xticks([]) + self.ax4.set_yticks([]) + self.ax4.set_title('Message') + for i in range(len(full_message)): + self.ax4.text(0.05, 0.95 - i * 0.05, + full_message[i], + verticalalignment='top', + horizontalalignment='left', + transform=self.ax4.transAxes, + color=color, + fontsize=15) + self.i_fig.canvas.draw() + + return + + def display_help_text(self): + """Shows static text on the top right subplot + + This will print static help text on the top right subplot of the + interactive window. + + Notes: + This is really hard to format and having a proper GUI should help + to have probably richer formatted text on the screen. + + """ + self.ax2.set_title('Help') + self.ax2.set_xticks([]) + self.ax2.set_yticks([]) + + self.ax2.text(1, 11, 'F1 or ?:', + fontsize=13) + + self.ax2.text(1.46, 11, 'Prints Help.', + fontsize=13) + + self.ax2.text(1, 10.5, 'F2 or f:', + fontsize=13) + + self.ax2.text(1.46, 10.5, 'Fit Wavelength Solution to points', + fontsize=13) + + self.ax2.text(1.46, 10, 'collected', + fontsize=13) + + self.ax2.text(1, 9.5, 'F3 or a:', + fontsize=13) + + self.ax2.text(1.46, 9.5, 'Find new lines, use when the solution', + fontsize=13) + + self.ax2.text(1.46, 9, 'is already decent.', + fontsize=13) + + self.ax2.text(1, 8.5, 'F4:', + fontsize=13) + + self.ax2.text(1.46, 8.5, 'Evaluate Solution', + fontsize=13) + + self.ax2.text(1, 8, 'F6 or l:', + fontsize=13) + + self.ax2.text(1.46, 8, 'Linearize Data', + fontsize=13) + + self.ax2.text(1, 7.5, 'd :', + fontsize=13) + + self.ax2.text(1.46, 7.5, 'Delete Closest Point', + fontsize=13) + + self.ax2.text(1, 7, 'Ctrl+d:', + fontsize=13) + + self.ax2.text(1.46, 7, 'Delete all recorded marks.', + fontsize=13) + + self.ax2.text(1, 6.5, 'Ctrl+t:', + fontsize=13) + + self.ax2.text(1.46, 6.5, 'Save as reference lamp template.', + fontsize=13) + + self.ax2.text(1, 6, 'Ctrl+z:', + fontsize=13) + + self.ax2.text(1.46, 6, 'Remove all automatic added points.', + fontsize=13) + + self.ax2.text(1.46, 5.5, 'Undo what F3 does.', + fontsize=13) + + self.ax2.text(1, 5, 'Middle Button or p:', + fontsize=13) + + self.ax2.text(1.46, 4.5, 'Finds and records line position', + fontsize=13) + + self.ax2.text(1, 4, 'Enter :', + fontsize=13) + + self.ax2.text(1.46, 4, 'Close Figure and apply wavelength', + fontsize=13) + + self.ax2.text(1.46, 3.5, 'solution', + fontsize=13) + + self.ax2.set_ylim((0, 12)) + self.ax2.set_xlim((0.95, 3.5)) diff --git a/goodman_pipeline/spectroscopy/redspec.py b/goodman_pipeline/spectroscopy/redspec.py index 3afeaf6e..f2b12af9 100755 --- a/goodman_pipeline/spectroscopy/redspec.py +++ b/goodman_pipeline/spectroscopy/redspec.py @@ -17,7 +17,6 @@ from .wavelength import WavelengthCalibration from ..core import (classify_spectroscopic_data, search_comp_group, - add_wcs_keys, identify_targets, trace_targets, extraction, @@ -172,6 +171,16 @@ def get_args(arguments=None): "discriminate usable targets. Default 3 times " "background level") + parser.add_argument('--interactive-wavelength', + action='store_true', + dest='interactive_wavelength', + help="Interactive Wavelength Calibration.") + + parser.add_argument('--linearize', + action='store_true', + dest='linearize', + help="Save as a linear wavelength solution, requires resampling.") + parser.add_argument('--save-plots', action='store_true', dest='save_plots', @@ -237,16 +246,16 @@ def get_args(arguments=None): return args -class MainApp(object): +class RedSpec(object): """Defines and initialize all important variables for processing the data - The MainApp class controls the way the night is organized for further + The RedSpec class controls the way the night is organized for further processing. It also sets the appropriate parameters that will allow for a smooth working in all the other modules. """ def __init__(self): - """Init method for MainApp class + """Init method for RedSpec class This method initializes the arguments for the class, if they are not provided it will get them. @@ -261,7 +270,7 @@ def __init__(self): self._pipeline_version = __version__ def __call__(self, args=None): - """Call method for the MainApp class + """Call method for the RedSpec class This method call the higher level functions in order to do the spectroscopic data reduction. @@ -295,7 +304,7 @@ def __call__(self, args=None): else: self.log.debug("Received non-empty data container.") - self.log.debug("Calling _run method for MainApp") + self.log.debug("Calling _run method for RedSpec") self._run(data_container=data_container, extraction_type=self.args.extraction_type, target_fit_model=self.args.target_fit_model, @@ -328,7 +337,9 @@ def _run(self, obj_groupby = object_group.groupby(['object']).size( ).reset_index().rename(columns={0: 'count'}) - self.log.info("Processing Science Target: " + self.log.debug("Here is where the process starts") + + self.log.info("Starting process for Science Target: " "{:s} with {:d} files." "".format(obj_groupby.iloc[0]['object'], obj_groupby.iloc[0]['count'])) @@ -342,7 +353,7 @@ def _run(self, comp_group = self.reference.check_comp_group(comp_group) if comp_group is None: - self.log.debug('Group does not have comparison lamps') + self.log.warning('Group does not have comparison lamps') if data_container.comp_groups is not None: self.log.debug('There are comparison lamp group ' 'candidates') @@ -380,16 +391,15 @@ def _run(self, file_path = os.path.join(full_path, spec_file) ccd = CCDData.read(file_path, unit=u.adu) ccd.header.set('GSP_PNAM', value=spec_file) - ccd = add_wcs_keys(ccd=ccd) # ccd.header['GSP_FNAM'] = spec_file if comp_group is not None and comp_ccd_list == []: for comp_file in comp_group.file.tolist(): + self.log.debug(f"Preparing comparison lamp file {comp_file} for processing.") comp_path = os.path.join(full_path, comp_file) comp_ccd = CCDData.read(comp_path, unit=u.adu) - comp_ccd = add_wcs_keys(ccd=comp_ccd) comp_ccd.header.set('GSP_PNAM', value=comp_file) comp_ccd_list.append(comp_ccd) @@ -399,8 +409,7 @@ def _run(self, 'Comp Group is None or comp list already exist') # identify - self.log.debug("Calling procedure for target " - "identification.") + self.log.info(f"Starting target identification for file {spec_file}") target_list = identify_targets( ccd=ccd, fit_model=target_fit_model, @@ -424,12 +433,20 @@ def _run(self, continue # if len(trace_list) > 0: + # Experimental interactive extraction + # if False: + # ie = InteractiveExtraction() + # # ie = InteractiveExtractionGUI() + # ie(ccd=ccd, lamps=comp_ccd_list, traces=trace_list) + extracted_target_and_lamps = [] for single_trace, single_profile, trace_info in trace_list: if single_profile.__class__.name == 'Gaussian1D': single_profile_center = single_profile.mean.value elif single_profile.__class__.name == 'Moffat1D': single_profile_center = single_profile.x_0.value + self.log.info(f"Processing traced profile fitted wit {single_profile.__class__.name} " + f"centered at {single_profile_center}") if len(trace_list) > 1: target_number = trace_list.index( @@ -454,6 +471,7 @@ def _run(self, # print(spec_file) # lamp extraction + self.log.debug(f"Starting lamps extraction for science file {spec_file}") all_lamps = [] if comp_ccd_list: for comp_lamp in comp_ccd_list: @@ -525,6 +543,7 @@ def _run(self, self.log.error('No target was identified in file' ' {:s}'.format(spec_file)) continue + self.log.debug(f"Starting Wavelength Calibration Process for File {spec_file}.") object_number = None for sci_target, comp_list in extracted_target_and_lamps: try: @@ -536,18 +555,22 @@ def _run(self, reference_data=self.args.reference_dir, object_number=object_number, output_prefix=self.args.output_prefix, + interactive_wavelength=self.args.interactive_wavelength, + linearize=self.args.linearize, plot_results=self.args.plot_results, save_plots=self.args.save_plots, plots=self.args.debug_with_plots) except NoMatchFound as no_match_error: self.log.error(no_match_error) - except NotImplemented as error: + except NotImplementedError as error: self.log.error(error) + self.log.debug(f"End of process for file {spec_file}.") + if __name__ == '__main__': # pragma: no cover - MAIN_APP = MainApp() + RED_SPEC = RedSpec() try: - MAIN_APP() + RED_SPEC() except KeyboardInterrupt: sys.exit(0) diff --git a/goodman_pipeline/spectroscopy/tests/__init__.py b/goodman_pipeline/spectroscopy/tests/__init__.py index 3be6bd29..c3961685 100644 --- a/goodman_pipeline/spectroscopy/tests/__init__.py +++ b/goodman_pipeline/spectroscopy/tests/__init__.py @@ -1 +1 @@ -from __future__ import absolute_import \ No newline at end of file +from __future__ import absolute_import diff --git a/goodman_pipeline/spectroscopy/tests/test_redspec.py b/goodman_pipeline/spectroscopy/tests/test_redspec.py index 534c846c..b9adeec6 100644 --- a/goodman_pipeline/spectroscopy/tests/test_redspec.py +++ b/goodman_pipeline/spectroscopy/tests/test_redspec.py @@ -1,11 +1,9 @@ from __future__ import absolute_import -import argparse import os -import shutil from unittest import TestCase -from ...spectroscopy.redspec import (get_args, MainApp) +from ...spectroscopy.redspec import (get_args, RedSpec) class TestArguments(TestCase): @@ -57,17 +55,6 @@ def test_absolute_path_exists(self): self.assertTrue(os.path.exists(args.source)) self.assertTrue(os.path.isdir(args.source)) - # If yes, carry on - - # If not, convert it to become absolute - - # Check if the source folder exists - - # If it exists, carry on - - # If source folder does not exists, print message and leave program - # error_message = 'Input folder "{} does not exists."' - def test_get_args_output_path_does_not_exist(self): arguments = ['--data-path', 'does_not_exist'] self.assertRaises(SystemExit, get_args, arguments) @@ -118,20 +105,20 @@ def test_get_args(): return args -class TestMainApp(TestCase): +class TestRedSpec(TestCase): def setUp(self): - self.main_app = MainApp() + self.red_spec = RedSpec() def test_instantiation_without_args(self): - self.assertIsInstance(self.main_app, MainApp) - self.assertIsNone(self.main_app.args) - self.assertIsNone(self.main_app.wavelength_solution_obj) - self.assertIsNone(self.main_app.wavelength_calibration) - self.assertIsNone(self.main_app.reference) + self.assertIsInstance(self.red_spec, RedSpec) + self.assertIsNone(self.red_spec.args) + self.assertIsNone(self.red_spec.wavelength_solution_obj) + self.assertIsNone(self.red_spec.wavelength_calibration) + self.assertIsNone(self.red_spec.reference) def test___call___no_args(self): - self.assertRaises(SystemExit, self.main_app) + self.assertRaises(SystemExit, self.red_spec) def test___call___with_valid_arguments(self): arguments = ['--data-path', './', @@ -140,7 +127,7 @@ def test___call___with_valid_arguments(self): '--output-prefix', 'w', '--extraction', 'fractional'] args = get_args(arguments=arguments) - self.assertRaises(SystemExit, self.main_app, args) + self.assertRaises(SystemExit, self.red_spec, args) if __name__ == '__main__': diff --git a/goodman_pipeline/spectroscopy/tests/test_wavelength.py b/goodman_pipeline/spectroscopy/tests/test_wavelength.py index e4b9ad3d..0e4eb6a9 100644 --- a/goodman_pipeline/spectroscopy/tests/test_wavelength.py +++ b/goodman_pipeline/spectroscopy/tests/test_wavelength.py @@ -1,20 +1,18 @@ from __future__ import absolute_import -import json import numpy as np import os import re -from astropy.convolution import convolve, Gaussian1DKernel, Box1DKernel from astropy.io import fits -from astropy.modeling import models, Model +from astropy.modeling import models from ccdproc import CCDData -from unittest import TestCase, skip +from unittest import TestCase from ..wavelength import (WavelengthCalibration) from ..redspec import get_args -from ...core import add_wcs_keys, write_fits -from ...core import ReferenceData, NoMatchFound +from ...core import add_linear_wcs_keys, write_fits +from ...core import NoMatchFound class WavelengthCalibrationTests(TestCase): @@ -35,7 +33,7 @@ def setUp(self): self.ccd = CCDData(data=np.random.random_sample(200), meta=fits.Header(), unit='adu') - self.ccd = add_wcs_keys(ccd=self.ccd) + self.ccd = add_linear_wcs_keys(ccd=self.ccd) self.ccd.header.set('SLIT', value='1.0_LONG_SLIT', comment="slit [arcsec]") @@ -158,7 +156,7 @@ def test___call___method_one_comparison_lamps(self): reference_data='goodman_pipeline/data/ref_comp', json_output=True) - self.assertEqual(json_output['error'], 'Unable to obtain wavelength solution') + self.assertEqual(json_output['error'], 'Unable to obtain wavelength solution.') self.assertEqual(json_output['warning'], '') self.assertEqual(json_output['wavelength_solution'], []) @@ -214,6 +212,3 @@ def test___call___method_two_solution(self): for _solution in json_output['wavelength_solution']: self.file_list.append(_solution['file_name']) self.file_list.append(_solution['reference_lamp']) - - - diff --git a/goodman_pipeline/spectroscopy/wavelength.py b/goodman_pipeline/spectroscopy/wavelength.py index 96cb1eb7..9e95479b 100644 --- a/goodman_pipeline/spectroscopy/wavelength.py +++ b/goodman_pipeline/spectroscopy/wavelength.py @@ -30,10 +30,13 @@ evaluate_wavelength_solution, get_lines_in_lamp, linearize_spectrum, + record_wavelength_solution_evaluation, write_fits) from ..core import (ReferenceData, NoMatchFound) +from ._interactive_wavelength import InteractiveWavelengthCalibration + log = logging.getLogger(__name__) @@ -73,6 +76,7 @@ def __init__(self): self.poly_order = 3 self.wcs = WCS() self.wsolution = None + self.linearize = False self.wcal_lamp_file = None self.sci_target_file = None self.n_points = None @@ -81,7 +85,7 @@ def __init__(self): self.cross_corr_tolerance = 5 self.reference_data_dir = None self.reference_data = None - self.calibration_lamp = '' + self.comparison_lamp_file_name = '' self.wcal_lamp_file = '' # Instrument configuration and spectral characteristics @@ -89,17 +93,19 @@ def __init__(self): self.parallel_binning = None def __call__(self, - ccd, - comp_list, - save_data_to, - reference_data, - object_number=None, - corr_tolerance=15, - output_prefix='w', - plot_results=False, - save_plots=False, - plots=False, - json_output=False): + ccd: CCDData, + comp_list: list, + save_data_to: str | os.PathLike, + reference_data: str, + object_number: int | None = None, + corr_tolerance: int = 15, + output_prefix: str = 'w', + interactive_wavelength: bool = False, + linearize: bool = False, + plot_results: bool = False, + save_plots: bool = False, + plots: bool = False, + json_output: bool = False): """Call method for the WavelengthSolution Class It takes extracted data and produces wavelength calibrated 1D FITS file. @@ -140,6 +146,8 @@ def __call__(self, assert isinstance(ccd, CCDData) assert isinstance(comp_list, list) + self.linearize = linearize + json_payload = {'wavelength_solution': [], 'warning': '', 'error': ''} @@ -155,77 +163,96 @@ def __call__(self, self.i_fig = None - log.info('Processing Science Target: ' - '{:s}'.format(ccd.header['OBJECT'])) + log.info(f"Starting Wavelength calibration of Science Target: {ccd.header['OBJECT']} " + f"File: {self.sci_target_file}.") if len(comp_list) == 0: log.warning("No comparison lamps were provided for file {}" "".format(self.sci_target_file)) log.error("Ending processing of {}".format(self.sci_target_file)) if json_output: - json_payload['error'] ='Unable to process without reference lamps' + json_payload['error'] = 'Unable to process without reference lamps' return json_payload else: return else: + log.debug(f"Science file {self.sci_target_file} has {len(comp_list)} comparison lamps.") wavelength_solutions = [] reference_lamp_names = [] for self.lamp in comp_list: - self.calibration_lamp = self.lamp.header['GSP_FNAM'] - log.info('Using reference lamp {}'.format(self.calibration_lamp)) + self.comparison_lamp_file_name = self.lamp.header['GSP_FNAM'] self.raw_pixel_axis = range(self.lamp.shape[0]) self.lamp_name = self.lamp.header['OBJECT'] - log.info('Processing Comparison Lamp: ' - '{:s}'.format(self.lamp_name)) + log.info(f"Using Comparison lamp {self.lamp_name} {self.comparison_lamp_file_name}") self.lines_center = get_lines_in_lamp( ccd=self.lamp, plots=plots) try: - self._automatic_wavelength_solution( - save_data_to=save_data_to, - corr_tolerance=self.cross_corr_tolerance) + if interactive_wavelength: + interactive_wavelength = InteractiveWavelengthCalibration() + interactive_wavelength(ccd=ccd, + comparison_lamp=self.lamp, + save_data_to=save_data_to, + reference_data=reference_data) + self.wsolution = interactive_wavelength.wavelength_solution + self.rms_error, self.n_points, self.n_rejections = interactive_wavelength.wavelength_solution_evaluation + + else: + self._automatic_wavelength_solution( + save_data_to=save_data_to, + corr_tolerance=self.cross_corr_tolerance) except NoMatchFound as message: - raise NoMatchFound(message) + log.error(message) + continue if self.wsolution is not None: - ccd.header.set('GSP_WRMS', value=self.rms_error) - ccd.header.set('GSP_WPOI', value=self.n_points) - ccd.header.set('GSP_WREJ', value=self.n_rejections) - - linear_x_axis, self.lamp.data = linearize_spectrum( - self.lamp.data, - wavelength_solution=self.wsolution) - - self.lamp = self.wcs.write_gsp_wcs(ccd=self.lamp, - model=self.wsolution) - - self.lamp = add_linear_wavelength_solution( - ccd=self.lamp, - x_axis=linear_x_axis, - reference_lamp=self.calibration_lamp) - - self.wcal_lamp_file = self._save_wavelength_calibrated( - ccd=self.lamp, - original_filename=self.calibration_lamp, - save_data_to=save_data_to, - output_prefix=output_prefix, - index=object_number, - lamp=True) + ccd = record_wavelength_solution_evaluation(ccd=ccd, + rms_error=self.rms_error, + n_points=self.n_points, + n_rejections=self.n_rejections) + self.lamp = record_wavelength_solution_evaluation(ccd=self.lamp, + rms_error=self.rms_error, + n_points=self.n_points, + n_rejections=self.n_rejections) + + if self.linearize: + + self.lamp.header.set('GSP_LINE', value='TRUE') + ccd.header.set('GSP_LINE', value='TRUE') + linear_x_axis, self.lamp.data = linearize_spectrum( + self.lamp.data, + wavelength_solution=self.wsolution) + + self.lamp = self.wcs.write_gsp_wcs(ccd=self.lamp, + model=self.wsolution) + + self.lamp = add_linear_wavelength_solution( + ccd=self.lamp, + x_axis=linear_x_axis, + reference_lamp=self.comparison_lamp_file_name) + + self.wcal_lamp_file = self._save_wavelength_calibrated( + ccd=self.lamp, + original_filename=self.comparison_lamp_file_name, + save_data_to=save_data_to, + output_prefix=output_prefix, + index=object_number, + lamp=True) + else: + pass wavelength_solutions.append(self.wsolution) reference_lamp_names.append(self.wcal_lamp_file) else: - log.error('It was not possible to get a wavelength ' - 'solution from lamp ' - '{:s} {:s}.'.format( - self.lamp.header['GSP_FNAM'], - self.lamp.header['OBJECT'])) + log.error(f"It was not possible to get a wavelength solution from lamp {self.lamp_name} " + f"{self.comparison_lamp_file_name}.") continue if len(wavelength_solutions) > 1: + log.warning(f"Multiple ({len(wavelength_solutions)}) wavelength solutions found.") warning_message = str("The current version of the pipeline " "does not combine multiple solution " "instead it saves a single version of " @@ -281,9 +308,9 @@ def __call__(self, return json_payload else: - log.error("No wavelength solution.") + log.error("Unable to obtain wavelength solution.") if json_output: - json_payload['error'] = "Unable to obtain wavelength solution" + json_payload['error'] = "Unable to obtain wavelength solution." return json_payload def _automatic_wavelength_solution(self, @@ -332,10 +359,11 @@ def _automatic_wavelength_solution(self, reference_lamp_ccd = self.reference_data.get_reference_lamp( header=self.lamp.header) - log.debug('Found reference lamp: ' - '{:s}'.format(reference_lamp_ccd.header['GSP_FNAM'])) - except NoMatchFound as error: - raise NoMatchFound(error) + if reference_lamp_ccd is not None: + log.debug(f"Found reference lamp {reference_lamp_ccd.header['GSP_FNAM']} " + f"for {self.sci_target_file}'s lamp {self.lamp.header['OBJECT']}: ") + else: + raise NoMatchFound(f"Unable to find a suitable reference lamp for {self.lamp.header['OBJECT']} {self.lamp.header['GSP_FNAM']}.") except NotImplementedError as error: raise NotImplementedError(error) @@ -377,7 +405,8 @@ def _automatic_wavelength_solution(self, reference=reference_lamp_ccd.data, compared=self.lamp.data, slit_size=slit_size, - serial_binning=self.serial_binning) + serial_binning=self.serial_binning, + selection_bias='center') log.debug(f"Found global cross-correlation value of: {global_cross_corr}") if plots: @@ -491,7 +520,7 @@ def _automatic_wavelength_solution(self, if self.wsolution is None: log.error('Failed to find wavelength solution using reference ' - 'file: {:s}'.format(self.calibration_lamp)) + 'file: {:s}'.format(self.comparison_lamp_file_name)) return None # finding differences in order to improve the wavelength solution @@ -662,14 +691,18 @@ def _save_science_data(self, """ ccd = ccd.copy() - linear_x_axis, ccd.data = linearize_spectrum( - data=ccd.data, - wavelength_solution=wavelength_solution) - - ccd = add_linear_wavelength_solution( - ccd=ccd, - x_axis=linear_x_axis, - reference_lamp=self.calibration_lamp) + resample = False + if resample: + linear_x_axis, ccd.data = linearize_spectrum( + data=ccd.data, + wavelength_solution=wavelength_solution) + + ccd = add_linear_wavelength_solution( + ccd=ccd, + x_axis=linear_x_axis, + reference_lamp=self.wcal_lamp_file) + else: + wavelength_axis = wavelength_solution(range(ccd.data.shape[0])) save_file_name = self._save_wavelength_calibrated( ccd=ccd, @@ -726,8 +759,7 @@ def _save_science_data(self, plot_path = os.path.join(plots_dir, plot_name) # print(plot_path) plt.savefig(plot_path, dpi=300) - log.info('Saved plot as {:s} file ' - 'DPI=300'.format(plot_name)) + log.info('Saved plot as {:s} file DPI=300'.format(plot_name)) if plots or plot_results: # pragma: no cover manager = plt.get_current_fig_manager() @@ -758,8 +790,7 @@ def _save_wavelength_calibrated(self, f_end = '_ws_{:d}.fits'.format(index) file_full_path = os.path.join(save_data_to, - output_prefix + - original_filename.replace('.fits', f_end)) + output_prefix + original_filename.replace('.fits', f_end)) if lamp: log.info('Wavelength-calibrated {:s} file saved to: ' diff --git a/goodman_pipeline/wcs/__init__.py b/goodman_pipeline/wcs/__init__.py index ce7a0599..1b72ac48 100644 --- a/goodman_pipeline/wcs/__init__.py +++ b/goodman_pipeline/wcs/__init__.py @@ -1 +1 @@ -from .wcs import WCS \ No newline at end of file +from .wcs import WCS diff --git a/goodman_pipeline/wcs/tests/test_functional.py b/goodman_pipeline/wcs/tests/test_functional.py index 706b0212..29923bc8 100644 --- a/goodman_pipeline/wcs/tests/test_functional.py +++ b/goodman_pipeline/wcs/tests/test_functional.py @@ -1,17 +1,18 @@ from __future__ import absolute_import -from unittest import TestCase, skip -from ..wcs import WCS +from unittest import TestCase + import numpy as np import os import re import sys -from astropy.io import fits -from astropy.modeling import (models, fitting, Model) -import matplotlib.pyplot as plt +from astropy.io import fits +from astropy.modeling import (models, Model) from ccdproc import CCDData +from ..wcs import WCS + class TestWCSBase(TestCase): @@ -54,7 +55,7 @@ def test_fit_chebyshev(self): self.assertEqual(model.degree, ccd.header['GSP_ORDR']) for i in range(model.degree + 1): self.assertAlmostEqual(model.__getattribute__('c{:d}'.format(i)).value, - ccd.header['GSP_C{:03d}'.format(i)]) + ccd.header['GSP_C{:03d}'.format(i)]) def test_fit_linear(self): test_file = os.path.join(self.data_path, @@ -127,7 +128,7 @@ def test_read__non_linear_chebyshev(self): ccd = CCDData.read(test_file, unit='adu') - result = self.wcs.read(ccd=ccd) + self.wcs.read(ccd=ccd) self.assertIsInstance(self.wcs.model, Model) self.assertEqual(self.wcs.model.__class__.__name__, 'Chebyshev1D') @@ -138,7 +139,7 @@ def test_read__non_linear_legendre(self): ccd = CCDData.read(test_file, unit='adu') - result = self.wcs.read(ccd=ccd) + self.wcs.read(ccd=ccd) self.assertIsInstance(self.wcs.model, Model) self.assertEqual(self.wcs.model.__class__.__name__, 'Legendre1D') @@ -193,8 +194,8 @@ def test_write_gsp_wcs(self): self.assertIsInstance(model, Model) blank_ccd = CCDData(data=np.ones(ccd.data.shape), - meta=fits.Header(), - unit='adu') + meta=fits.Header(), + unit='adu') blank_ccd.header.set('GSP_WREJ', value=None, comment='empty') new_ccd = self.wcs.write_gsp_wcs(ccd=blank_ccd, model=model) @@ -204,7 +205,7 @@ def test_write_gsp_wcs(self): self.assertEqual(new_ccd.header['GSP_NPIX'], ccd.header['GSP_NPIX']) for i in range(model.degree + 1): self.assertAlmostEqual(new_ccd.header['GSP_C{:03d}'.format(i)], - ccd.header['GSP_C{:03d}'.format(i)]) + ccd.header['GSP_C{:03d}'.format(i)]) def test_read_gsp_wcs(self): test_file = os.path.join(self.data_path, diff --git a/goodman_pipeline/wcs/tests/test_unittest.py b/goodman_pipeline/wcs/tests/test_unittest.py index 746363bd..d6130da6 100644 --- a/goodman_pipeline/wcs/tests/test_unittest.py +++ b/goodman_pipeline/wcs/tests/test_unittest.py @@ -1,19 +1,13 @@ from __future__ import absolute_import -from unittest import TestCase, skip +from unittest import TestCase from ..wcs import WCS -import numpy as np -import os -import re -import sys -from astropy.io import fits -from astropy.modeling import (models, fitting, Model) -from ccdproc import CCDData +from astropy.modeling import Model class TestWCS(TestCase): - + def setUp(self): self.wcs = WCS() @@ -66,7 +60,7 @@ def test_pm_fitter_undefined_model_and_fitter(self): # self.wcs._fitter(physical=pixel, wavelength=angstrom) def test_pm_fitter_not_enough_points(self): - pixel = [1,2] + pixel = [1, 2] angstrom = [8000, 8001] self.wcs.model_name = 'chebyshev' @@ -76,7 +70,6 @@ def test_pm_fitter_not_enough_points(self): result = self.wcs._fitter(physical=pixel, wavelength=angstrom) self.assertIsNone(result) - def test_pm_set_math_model__none(self): self.wcs.wcs_dict['dtype'] = -1 self.assertRaises(NotImplementedError, self.wcs._set_math_model) @@ -113,7 +106,7 @@ def test_pm_set_math_model__wrong_ftype(self): self.assertRaisesRegex(SyntaxError, 'ftype {:d} is not defined in the ' 'standard'.format(self.wcs.wcs_dict['ftype'])) - + def test_pm_none(self): self.assertRaises(NotImplementedError, self.wcs._none) @@ -155,4 +148,3 @@ def test_pm_non_linear_lspline(self): def test_pm_non_linear_cspline(self): self.assertRaises(NotImplementedError, self.wcs._non_linear_cspline) - diff --git a/goodman_pipeline/wcs/wcs.py b/goodman_pipeline/wcs/wcs.py index 21c0cb27..a6431c47 100644 --- a/goodman_pipeline/wcs/wcs.py +++ b/goodman_pipeline/wcs/wcs.py @@ -252,7 +252,7 @@ def _fitter(self, physical, wavelength): # wavelength solution reader private methods. def _read_non_linear(self, dimension): - """Non linear solutions reader + """Non-linear solutions reader Notes: Not all kind of non-linear solutions are implemented. Apparently is diff --git a/requirements.txt b/requirements.txt index 43dc3f23..40d1f507 100644 --- a/requirements.txt +++ b/requirements.txt @@ -2,6 +2,7 @@ astropy astroplan astroscrappy ccdproc +coloredlogs coveralls cython recommonmark @@ -9,6 +10,7 @@ matplotlib mock numpy pandas +setuptools scipy sphinx sphinxcontrib.napoleon diff --git a/setup.py b/setup.py index caa5c8ed..61afabe1 100644 --- a/setup.py +++ b/setup.py @@ -155,6 +155,4 @@ def create_version_py(packagename, version, source_dir='.'): scripts=['goodman_pipeline/scripts/redccd', 'goodman_pipeline/scripts/redspec', ], - ) - - +)