From e6a36d44f5733dc72c83c38bcf1e1b5ec7a564cf Mon Sep 17 00:00:00 2001 From: Alfredo Mejia-Narvaez Date: Thu, 16 Jan 2025 10:49:15 -0300 Subject: [PATCH 1/3] improving routine for extraction validation --- python/lvmdrp/functions/imageMethod.py | 90 +++++++++++++++++--------- 1 file changed, 60 insertions(+), 30 deletions(-) diff --git a/python/lvmdrp/functions/imageMethod.py b/python/lvmdrp/functions/imageMethod.py index 8c000709..241d8376 100644 --- a/python/lvmdrp/functions/imageMethod.py +++ b/python/lvmdrp/functions/imageMethod.py @@ -2975,52 +2975,82 @@ def reprojectRSS_drp( rep.writeto(f"{out_path}/{out_name}_2d.fits", overwrite=True) -def testres_drp(image, trace, fwhm, flux): - """ - Historic task used for debugging of the the extraction routine... +def validate_extraction(in_image, in_cent, in_width, in_rss): + """Evaluates the extracted flux into the original 2D pixel grid + + This routine will evaluate the extracted flux in the original + 2D grid and compare the resulting 2D model against the original + 2D image. Three images are stored as outputs: + + - The residual 2D image: model - data + - The ratio 2D image: model / data + - The 2D model + + Parameters + ---------- + in_image : str + path to the original 2D image of the extracted flux + in_cent : str + path to the fiber centroids trace + in_width : str + path to the fiber width (Fn_iWHM) trace + flux : str + path to the extracted flux in RSS format """ + log.info(f"loading 2D image {in_image}") img = Image() - # t1 = time.time() - img.loadFitsData(image, extension_data=0) - trace_mask = TraceMask() - trace_mask.loadFitsData(trace, extension_data=0) - trace_fwhm = TraceMask() - # trace_fwhm.setData(data=numpy.ones(trace_mask._data.shape)*2.5) - trace_fwhm.loadFitsData(fwhm, extension_data=0) + img._data = numpy.nan_to_num(img._data) + img.loadFitsData(in_image, extension_data=0) + log.info(f"loading fiber parameters in {in_cent} and {in_width}") + cent = TraceMask.from_file(in_cent, extension_data=0) + width = TraceMask.from_file(in_width, extension_data=0) + width._data /= 2.354 + + log.info(f"loading extracted flux in {in_rss}") trace_flux = TraceMask() - trace_flux.loadFitsData(flux, extension_data=0) + trace_flux.loadFitsData(in_rss, extension_data=0) + trace_flux._data = numpy.nan_to_num(trace_flux._data) + + ypix_cor = trace_flux._slitmap[["spectrographid"] == int(img._header["CCD"][1])]["ypix_z"] + ypix_ori = img._slitmap[["spectrographid"] == int(img._header["CCD"][1])]["ypix_z"] + thermal_shift = ypix_cor - ypix_ori + log.info(f"fiber thermal shift in slitmap: {thermal_shift[0]:.4f}") + cent._data += thermal_shift[:, None] + + log.info(f"evaluating extracted flux into 2D pixel grid for {img._dim[1]} columns") x = numpy.arange(img._dim[0]) out = numpy.zeros(img._dim) fact = numpy.sqrt(2.0 * numpy.pi) for i in range(img._dim[1]): # print i - A = ( - 1.0 - * numpy.exp( - -0.5 - * ( - (x[:, numpy.newaxis] - trace_mask._data[:, i][numpy.newaxis, :]) - / abs(trace_fwhm._data[:, i][numpy.newaxis, :] / 2.354) - ) - ** 2 - ) - / (fact * abs(trace_fwhm._data[:, i][numpy.newaxis, :] / 2.354)) - ) + A = (numpy.exp(-0.5 * ((x[:, None] - cent._data[:, i][None, :]) / abs(width._data[:, i][None, :])) ** 2) / (fact * abs(width._data[:, i][None, :]))) + # print(numpy.isnan(A).sum(), numpy.isnan(trace_flux._data[:, i]).sum()) spec = numpy.dot(A, trace_flux._data[:, i]) out[:, i] = spec if i == 1000: - plt.plot(spec, "-r") - plt.plot(img._data[:, i], "ok") + plt.figure() + plt.step(x, img._datin_a[:, i], color="k", lw=1, where="mid") + plt.step(x, spec, color="r", lw=1, where="mid") plt.show() + out_path = os.path.dirname(in_image) + out_name = os.path.basename(in_image).split(".fits")[0] + out_residual = os.path.join(out_path, f"{out_name}_residual.fits") + out_2dimage = os.path.join(out_path, f"{out_name}_2dimage.fits") + out_ratio = os.path.join(out_path, f"{out_name}_ratio.fits") + + log.info(f"writing residual to {out_residual}") hdu = pyfits.PrimaryHDU(img._data - out) - hdu.writeto("res.fits", overwrite=True) - hdu = pyfits.PrimaryHDU(out) - hdu.writeto("fit.fits", overwrite=True) + hdu.writeto(out_residual, overwrite=True) - hdu = pyfits.PrimaryHDU((img._data - out) / img._data) - hdu.writeto("res_rel.fits", overwrite=True) + log.info(f"writing ratio to {out_ratio}") + hdu = pyfits.PrimaryHDU(out / img._data) + hdu.writeto(out_ratio, overwrite=True) + + log.info(f"writing 2D model {out_2dimage}") + hdu = pyfits.PrimaryHDU(out) + hdu.writeto(out_2dimage, overwrite=True) # TODO: for arcs take short exposures for bright lines & long exposures for faint lines From 2870946a943b4c45f261968d1f7285bddda23135 Mon Sep 17 00:00:00 2001 From: Alfredo Mejia-Narvaez Date: Thu, 16 Jan 2025 12:09:21 -0300 Subject: [PATCH 2/3] fixed busg in validate_extraction routine & improved QA plots --- python/lvmdrp/functions/imageMethod.py | 48 ++++++++++++++------------ 1 file changed, 25 insertions(+), 23 deletions(-) diff --git a/python/lvmdrp/functions/imageMethod.py b/python/lvmdrp/functions/imageMethod.py index 241d8376..8eb1d1b2 100644 --- a/python/lvmdrp/functions/imageMethod.py +++ b/python/lvmdrp/functions/imageMethod.py @@ -2975,7 +2975,7 @@ def reprojectRSS_drp( rep.writeto(f"{out_path}/{out_name}_2d.fits", overwrite=True) -def validate_extraction(in_image, in_cent, in_width, in_rss): +def validate_extraction(in_image, in_cent, in_width, in_rss, plot_columns=[1000, 2000, 3000], display_plots=False): """Evaluates the extracted flux into the original 2D pixel grid This routine will evaluate the extracted flux in the original @@ -2989,56 +2989,58 @@ def validate_extraction(in_image, in_cent, in_width, in_rss): Parameters ---------- in_image : str - path to the original 2D image of the extracted flux + Path to the original 2D image of the extracted flux in_cent : str - path to the fiber centroids trace + Path to the fiber centroids trace in_width : str - path to the fiber width (Fn_iWHM) trace - flux : str - path to the extracted flux in RSS format + Path to the fiber width (FWHM) trace + in_rss : str + Path to the extracted flux in RSS format + plot_columns : array_like, optional + columns to show in plot, by default [1000, 2000, 3000] + display_plots : bool, optional + whether to display plots to screen or not, by dafult False """ log.info(f"loading 2D image {in_image}") img = Image() img._data = numpy.nan_to_num(img._data) - img.loadFitsData(in_image, extension_data=0) + img.loadFitsData(in_image) log.info(f"loading fiber parameters in {in_cent} and {in_width}") - cent = TraceMask.from_file(in_cent, extension_data=0) - width = TraceMask.from_file(in_width, extension_data=0) + cent = TraceMask.from_file(in_cent) + width = TraceMask.from_file(in_width) width._data /= 2.354 log.info(f"loading extracted flux in {in_rss}") - trace_flux = TraceMask() - trace_flux.loadFitsData(in_rss, extension_data=0) - trace_flux._data = numpy.nan_to_num(trace_flux._data) + rss = RSS.from_file(in_rss) + rss._data = numpy.nan_to_num(rss._data) - ypix_cor = trace_flux._slitmap[["spectrographid"] == int(img._header["CCD"][1])]["ypix_z"] + ypix_cor = rss._slitmap[["spectrographid"] == int(img._header["CCD"][1])]["ypix_z"] ypix_ori = img._slitmap[["spectrographid"] == int(img._header["CCD"][1])]["ypix_z"] thermal_shift = ypix_cor - ypix_ori - log.info(f"fiber thermal shift in slitmap: {thermal_shift[0]:.4f}") - cent._data += thermal_shift[:, None] + log.info(f"fiber thermal shift in slitmap: {thermal_shift:.4f}") + cent._data += thermal_shift log.info(f"evaluating extracted flux into 2D pixel grid for {img._dim[1]} columns") x = numpy.arange(img._dim[0]) out = numpy.zeros(img._dim) fact = numpy.sqrt(2.0 * numpy.pi) + + fig, axs = create_subplots(to_display=display_plots, nrows=len(plot_columns), ncols=1, figsize=(15,5), sharex=True, layout="constrained") for i in range(img._dim[1]): - # print i A = (numpy.exp(-0.5 * ((x[:, None] - cent._data[:, i][None, :]) / abs(width._data[:, i][None, :])) ** 2) / (fact * abs(width._data[:, i][None, :]))) - # print(numpy.isnan(A).sum(), numpy.isnan(trace_flux._data[:, i]).sum()) - spec = numpy.dot(A, trace_flux._data[:, i]) + spec = numpy.dot(A, rss._data[:, i]) out[:, i] = spec - if i == 1000: - plt.figure() - plt.step(x, img._datin_a[:, i], color="k", lw=1, where="mid") - plt.step(x, spec, color="r", lw=1, where="mid") - plt.show() + if i in plot_columns: + axs[plot_columns.index(i)].step(x, img._data[:, i], color="k", lw=1, where="mid") + axs[plot_columns.index(i)].step(x, spec, color="r", lw=1, where="mid") out_path = os.path.dirname(in_image) out_name = os.path.basename(in_image).split(".fits")[0] out_residual = os.path.join(out_path, f"{out_name}_residual.fits") out_2dimage = os.path.join(out_path, f"{out_name}_2dimage.fits") out_ratio = os.path.join(out_path, f"{out_name}_ratio.fits") + save_fig(fig, product_path=out_2dimage, to_display=display_plots, figure_path="qa", label="2D_extracted_model") log.info(f"writing residual to {out_residual}") hdu = pyfits.PrimaryHDU(img._data - out) From d502cf447470f55749c0a0863329a68b026d384e Mon Sep 17 00:00:00 2001 From: Alfredo Mejia-Narvaez Date: Fri, 24 Jan 2025 15:17:43 -0300 Subject: [PATCH 3/3] fixing Halpha/Hbeta/Hgamma in sensitivities --- python/lvmdrp/core/fluxcal.py | 4 ++-- python/lvmdrp/functions/fluxCalMethod.py | 22 ++++++++++++++++++++-- 2 files changed, 22 insertions(+), 4 deletions(-) diff --git a/python/lvmdrp/core/fluxcal.py b/python/lvmdrp/core/fluxcal.py index a49ff8b4..5d1e2bb3 100644 --- a/python/lvmdrp/core/fluxcal.py +++ b/python/lvmdrp/core/fluxcal.py @@ -172,9 +172,9 @@ def butter_lowpass_filter(data, cutoff_freq, nyq_freq, order=4): return y -def filter_channel(w, f, k=3): +def filter_channel(w, f, k=3, cutoff_freq=0.01): c = np.where(np.isfinite(f)) - s = butter_lowpass_filter(f[c], 0.01, 2) + s = butter_lowpass_filter(f[c], cutoff_freq, 2) res = s - f[c] # plt.plot(w[c], f[c], 'k.') # plt.plot(w[c], s, 'b-') diff --git a/python/lvmdrp/functions/fluxCalMethod.py b/python/lvmdrp/functions/fluxCalMethod.py index e1036a0c..bbbc6c20 100644 --- a/python/lvmdrp/functions/fluxCalMethod.py +++ b/python/lvmdrp/functions/fluxCalMethod.py @@ -222,13 +222,15 @@ def apply_fluxcal(in_rss: str, out_fframe: str, method: str = 'STD', display_plo return fframe -def standard_sensitivity(stds, rss, GAIA_CACHE_DIR, ext, res, plot=False, width=3): +def standard_sensitivity(stds, rss, GAIA_CACHE_DIR, ext, res, plot=False, width=3, cutoff_freq=0.01): # load the sky masks channel = rss._header['CCD'] w = rss._wave m = get_sky_mask_uves(w, width=width) m2 = None + # if channel == "r": + # m2 = ~(((w>6563-80)&(w<6563+80))|((w>6850)&(w<6950))) if channel == "z": m2 = get_z_continuum_mask(w) @@ -262,6 +264,16 @@ def standard_sensitivity(stds, rss, GAIA_CACHE_DIR, ext, res, plot=False, width= spec = (rss._data[fibidx[0],:] - master_sky._data[fibidx[0],:])/exptime #TODO: check exptime? + # remove nans + # spec[(((w>6563-80)&(w<6563+80))|((w>6850)&(w<6950)))] = np.nan + # plt.figure() + # plt.plot(w, spec) + + spec = fluxcal.interpolate_mask(w, spec, np.isnan(spec), fill_value="extrapolate") + + # print(">>>>>>>>>>", w[(((w>6563-80)&(w<6563+80))|((w>6850)&(w<6950)))]) + + # interpolate over bright sky lines spec = fluxcal.interpolate_mask(w, spec, m, fill_value="extrapolate") if channel == "z": @@ -277,7 +289,13 @@ def standard_sensitivity(stds, rss, GAIA_CACHE_DIR, ext, res, plot=False, width= # divide to find sensitivity and smooth sens = stdflux / spec - wgood, sgood = fluxcal.filter_channel(w, sens, 2) + if channel == "r": + sens[(((w>6563-80)&(w<6563+80))|((w>6850)&(w<6950)))] = np.nan + if channel == "b": + sens[(((w>4860-50)&(w<4860+50))|((w>4340-50)&(w<4340+30)))] = np.nan + wgood, sgood = fluxcal.filter_channel(w, sens, 3, cutoff_freq=cutoff_freq) + # wgood = w[np.isfinite(sens)] + # sgood = sens[np.isfinite(sens)] s = interpolate.make_smoothing_spline(wgood, sgood, lam=1e4) res[f"STD{nn}SEN"] = s(w).astype(np.float32) * u.Unit("erg / (ct cm2)")