diff --git a/python/lvmdrp/core/image.py b/python/lvmdrp/core/image.py index 540610cc..159d5e65 100644 --- a/python/lvmdrp/core/image.py +++ b/python/lvmdrp/core/image.py @@ -23,7 +23,7 @@ from lvmdrp.core.apertures import Apertures from lvmdrp.core.header import Header from lvmdrp.core.tracemask import TraceMask -from lvmdrp.core.spectrum1d import Spectrum1D, _cross_match_float, _cross_match, _spec_from_lines +from lvmdrp.core.spectrum1d import Spectrum1D, _normalize_peaks, _cross_match_float, _cross_match, _spec_from_lines from cextern.fast_median.fast_median import fast_median_filter_2d @@ -714,7 +714,7 @@ def add_header_comment(self, comstr): ''' self._header.append(('COMMENT', comstr), bottom=True) - def measure_fiber_shifts(self, ref_image, columns=[500, 1000, 1500, 2000, 2500, 3000], column_width=25, shift_range=[-5,5], axs=None): + def measure_fiber_shifts(self, ref_image, trace_cent, columns=[500, 1000, 1500, 2000, 2500, 3000], column_width=25, shift_range=[-5,5], axs=None): '''Measure the (thermal, flexure, ...) shift between the fiber (traces) in 2 detrended images in the y (cross dispersion) direction. Uses cross-correlations between (medians of a number of) columns to determine @@ -742,20 +742,42 @@ def measure_fiber_shifts(self, ref_image, columns=[500, 1000, 1500, 2000, 2500, elif isinstance(ref_image, numpy.ndarray): ref_data = ref_image + # unpack axes + axs_cc, axs_fb = axs + shifts = numpy.zeros(len(columns)) + select_blocks = [9] for j,c in enumerate(columns): s1 = bn.nanmedian(ref_data[50:-50,c-column_width:c+column_width], axis=1) s2 = bn.nanmedian(self._data[50:-50,c-column_width:c+column_width], axis=1) snr = numpy.sqrt(bn.nanmedian(self._data[50:-50,c-column_width:c+column_width], axis=1)) + median_snr = bn.nanmedian(snr) min_snr = 5.0 - if bn.nanmedian(snr) > min_snr: - _, shifts[j], _ = _cross_match_float(s1, s2, numpy.array([1.0]), shift_range, gauss_window=[-3,3], min_peak_dist=5.0, ax=axs[j]) - else: - comstr = f"low SNR (<={min_snr}) for thermal shift at column {c}: {bn.nanmedian(snr):.4f}, assuming = 0.0" + if median_snr <= min_snr: + comstr = f"low SNR (<={min_snr}) for thermal shift at column {c}: {median_snr:.4f}, assuming = 0.0" log.warning(comstr) self.add_header_comment(comstr) shifts[j] = 0.0 + continue + + _, shifts[j], _ = _cross_match_float(s1, s2, numpy.array([1.0]), shift_range, gauss_window=[-3,3], min_peak_dist=5.0, ax=axs_cc[j]) + + blocks_pos = numpy.asarray(numpy.split(trace_cent._data[:, c], 18))[select_blocks] + blocks_bounds = [(int(bpos.min())-5, int(bpos.max())+5) for bpos in blocks_pos] + + for i, (bmin, bmax) in enumerate(blocks_bounds): + x = numpy.arange(bmax-bmin) + i*(bmax-bmin) + 10 + y_model = bn.nanmedian(ref_data[bmin:bmax, c-column_width:c+column_width], axis=1) + y_data = bn.nanmedian(self._data[bmin:bmax, c-column_width:c+column_width], axis=1) + y_model = _normalize_peaks(y_model, min_peak_dist=5.0) + y_data = _normalize_peaks(y_data, min_peak_dist=5.0) + axs_fb[j].step(x, y_data, color="0.2", lw=1.5, label="data" if i == 0 else None) + axs_fb[j].step(x, y_model, color="tab:blue", lw=1, label="model" if i == 0 else None) + axs_fb[j].step(x+shifts[j], numpy.interp(x+shifts[j], x, y_model), color="tab:red", lw=1, label="corr. model" if i == 0 else None) + axs_fb[j].set_title(f"measured shift {shifts[j]:.4f} pixel @ column {c} with SNR = {median_snr:.2f}") + axs_fb[j].set_ylim(-0.05, 1.3) + axs_fb[0].legend(loc=1, frameon=False, ncols=3) return shifts diff --git a/python/lvmdrp/functions/imageMethod.py b/python/lvmdrp/functions/imageMethod.py index 8faef37d..72ca359a 100644 --- a/python/lvmdrp/functions/imageMethod.py +++ b/python/lvmdrp/functions/imageMethod.py @@ -43,7 +43,7 @@ ) from lvmdrp.core.plot import plt, create_subplots, plot_detrend, plot_strips, plot_image_shift, plot_fiber_thermal_shift, save_fig from lvmdrp.core.rss import RSS -from lvmdrp.core.spectrum1d import Spectrum1D, _spec_from_lines, _normalize_peaks, _cross_match +from lvmdrp.core.spectrum1d import Spectrum1D, _spec_from_lines, _cross_match from lvmdrp.core.tracemask import TraceMask from lvmdrp.utils.hdrfix import apply_hdrfix from lvmdrp.utils.convert import dateobs_to_sjd, correct_sjd @@ -379,16 +379,13 @@ def _fix_fiber_thermal_shifts(image, trace_cent, trace_width=None, trace_amp=Non # generate the continuum model using the master traces only along the specific columns if fiber_model is None: fiber_model, _ = image.eval_fiber_model(trace_cent, trace_width, trace_amp=trace_amp, columns=columns, column_width=column_width) - fiber_model.writeFitsData("./test_model.fits") mjd = image._header["SMJD"] expnum = image._header["EXPOSURE"] camera = image._header["CCD"] - # unpack axes - axs_cc, axs_fb = axs # calculate thermal shifts - column_shifts = image.measure_fiber_shifts(fiber_model, columns=columns, column_width=column_width, shift_range=shift_range, axs=axs_cc) + column_shifts = image.measure_fiber_shifts(fiber_model, trace_cent, columns=columns, column_width=column_width, shift_range=shift_range, axs=axs) # shifts stats median_shift = bn.nanmedian(column_shifts, axis=0) std_shift = bn.nanstd(column_shifts, axis=0) @@ -404,42 +401,6 @@ def _fix_fiber_thermal_shifts(image, trace_cent, trace_width=None, trace_amp=Non trace_cent_fixed._coeffs[:, 0] += median_shift trace_cent_fixed.eval_coeffs() - select_blocks = [9] - for j, c in enumerate(columns): - blocks_pos = numpy.asarray(numpy.split(trace_cent._data[:, c], 18))[select_blocks] - blocks_bounds = [(int(bpos.min())-5, int(bpos.max())+5) for bpos in blocks_pos] - - for i, (bmin, bmax) in enumerate(blocks_bounds): - x = numpy.arange(bmax-bmin) + i*(bmax-bmin) + 10 - # y_models = fiber_model._data[bmin:bmax,c-column_width:c+column_width] - y_model = bn.nanmedian(fiber_model._data[bmin:bmax,c-column_width:c+column_width], axis=1) - y_data = bn.nanmedian(image._data[bmin:bmax, c-column_width:c+column_width], axis=1) - snr = numpy.sqrt(y_data.mean()) - y_model = _normalize_peaks(y_model, min_peak_dist=5.0) - y_data = _normalize_peaks(y_data, min_peak_dist=5.0) - # axs_fib[j].step(x, y_models * norm, color="0.7", lw=0.7, alpha=0.3) - axs_fb[j].step(x, y_data, color="0.2", lw=1.5, label="data" if i == 0 else None) - axs_fb[j].step(x, y_model, color="tab:blue", lw=1, label="model" if i == 0 else None) - axs_fb[j].step(x+column_shifts[j], numpy.interp(x+column_shifts[j], x, y_model), color="tab:red", lw=1, label="corr. model" if i == 0 else None) - axs_fb[j].set_title(f"measured shift {column_shifts[j]:.4f} pixel @ column {c} with SNR = {snr:.2f}") - axs_fb[j].set_ylim(-0.05, 1.3) - axs_fb[0].legend(loc=1, frameon=False, ncols=3) - - # deltas = TraceMask(data=numpy.zeros_like(trace_cent._data), mask=numpy.ones_like(trace_cent._data, dtype=bool)) - # deltas._data[:, columns] = column_shifts - # deltas._mask[:, columns] = False - # deltas.fit_polynomial(deg=4) - - # # fig, ax = create_subplots(to_display=True, figsize=(15,5)) - # # ax.plot(deltas._data[:, columns]-column_shifts, ".k") - - # trace_cent_fixed = copy(trace_cent) - # for ifiber in range(trace_cent._data.shape[0]): - # poly_trace = numpy.polynomial.Polynomial(trace_cent._coeffs[ifiber]) - # poly_deltas = numpy.polynomial.Polynomial(deltas._coeffs[ifiber]) - # trace_cent_fixed._coeffs[ifiber] = (poly_trace + poly_deltas).coef - # trace_cent_fixed.eval_coeffs() - return trace_cent_fixed, column_shifts, fiber_model