From 4ac47f8e3d95eafc31262b51269c2d02cc26a0dc Mon Sep 17 00:00:00 2001 From: Thibault Barnouin Date: Fri, 14 Mar 2025 14:30:30 +0100 Subject: [PATCH 01/33] Save the raw total flux image as PrimaryHDU --- package/FOC_reduction.py | 34 ++++++++++++++++++++++----------- package/lib/fits.py | 41 ++++++++++++++++++++++++++++++++++------ 2 files changed, 58 insertions(+), 17 deletions(-) diff --git a/package/FOC_reduction.py b/package/FOC_reduction.py index eab1da7..eb73dc9 100755 --- a/package/FOC_reduction.py +++ b/package/FOC_reduction.py @@ -40,13 +40,13 @@ def main(target=None, proposal_id=None, infiles=None, output_dir="./data", crop= display_crop = False # Background estimation - error_sub_type = "scott" # sqrt, sturges, rice, scott, freedman-diaconis (default) or shape (example (51, 51)) - subtract_error = 2.0 + error_sub_type = "scott" # sqrt, sturges, rice, freedman-diaconis, scott (default) or shape (example (51, 51)) + subtract_error = 3.0 display_bkg = True # Data binning - pxsize = 0.05 - pxscale = "arcsec" # pixel, arcsec or full + pxsize = 40 + pxscale = "px" # pixel, arcsec or full rebin_operation = "sum" # sum or average # Alignement @@ -59,15 +59,15 @@ def main(target=None, proposal_id=None, infiles=None, output_dir="./data", crop= # Smoothing smoothing_function = "combine" # gaussian_after, weighted_gaussian_after, gaussian, weighted_gaussian or combine - smoothing_FWHM = 0.075 # If None, no smoothing is done - smoothing_scale = "arcsec" # pixel or arcsec + smoothing_FWHM = 1.5 # If None, no smoothing is done + smoothing_scale = "px" # pixel or arcsec # Rotation rotate_North = True # Polarization map output P_cut = 5 # if >=1.0 cut on the signal-to-noise else cut on the confidence level in Q, U - SNRi_cut = 5.0 # I measurments with SNR>30, which implies an uncertainty in P of 4.7%. + SNRi_cut = 1.0 # I measurments with SNR>30, which implies an uncertainty in P of 4.7%. flux_lim = None # lowest and highest flux displayed on plot, defaults to bkg and maximum in cut if None scale_vec = 3 step_vec = 1 # plot all vectors in the array. if step_vec = 2, then every other vector will be plotted if step_vec = 0 then all vectors are displayed at full length @@ -182,6 +182,14 @@ def main(target=None, proposal_id=None, infiles=None, output_dir="./data", crop= norm=LogNorm(vmin=data_array[data_array > 0.0].min() * headers[0]["photflam"], vmax=data_array[data_array > 0.0].max() * headers[0]["photflam"]), ) + flux_data, flux_error, flux_mask, flux_head = ( + deepcopy(data_array.sum(axis=0)), + deepcopy(np.sqrt(np.sum(error_array**2, axis=0))), + deepcopy(data_mask), + deepcopy(headers[0]), + ) + flux_head["EXPTIME"] = np.sum([head["EXPTIME"] for head in headers]) + # Rebin data to desired pixel size. if (pxsize is not None) and not (pxsize == 1 and pxscale.lower() in ["px", "pixel", "pixels"]): data_array, error_array, headers, Dxy, data_mask = proj_red.rebin_array( @@ -233,6 +241,8 @@ def main(target=None, proposal_id=None, infiles=None, output_dir="./data", crop= I_bkg, Q_bkg, U_bkg, S_cov_bkg, data_mask_bkg, header_bkg = proj_red.rotate_Stokes( I_bkg, Q_bkg, U_bkg, S_cov_bkg, np.array(True).reshape(1, 1), header_bkg, SNRi_cut=None ) + flux_data, flux_error, flux_mask, flux_head = proj_red.rotate_data(np.array([flux_data]), np.array([flux_error]), flux_mask, [flux_head]) + flux_data, flux_error, flux_head = flux_data[0], flux_error[0], flux_head[0] # Compute polarimetric parameters (polarization degree and angle). P, debiased_P, s_P, s_P_P, PA, s_PA, s_PA_P = proj_red.compute_pol(I_stokes, Q_stokes, U_stokes, Stokes_cov, header_stokes) @@ -258,8 +268,10 @@ def main(target=None, proposal_id=None, infiles=None, output_dir="./data", crop= figname, data_folder=data_folder, return_hdul=True, + flux_data=flux_data, + flux_head=flux_head, ) - outfiles.append("/".join([data_folder, Stokes_hdul[0].header["FILENAME"] + ".fits"])) + outfiles.append("/".join([data_folder, Stokes_hdul["I_STOKES"].header["FILENAME"] + ".fits"])) # Step 5: # crop to desired region of interest (roi) @@ -269,15 +281,15 @@ def main(target=None, proposal_id=None, infiles=None, output_dir="./data", crop= stokescrop.crop() stokescrop.write_to("/".join([data_folder, figname + ".fits"])) Stokes_hdul, header_stokes = stokescrop.hdul_crop, stokescrop.hdul_crop[0].header - outfiles.append("/".join([data_folder, Stokes_hdul[0].header["FILENAME"] + ".fits"])) + outfiles.append("/".join([data_folder, Stokes_hdul["I_STOKES"].header["FILENAME"] + ".fits"])) data_mask = Stokes_hdul["data_mask"].data.astype(bool) print( "F_int({0:.0f} Angs) = ({1} ± {2})e{3} ergs.cm^-2.s^-1.Angs^-1".format( header_stokes["PHOTPLAM"], *sci_not( - Stokes_hdul[0].data[data_mask].sum() * header_stokes["PHOTFLAM"], - np.sqrt(Stokes_hdul[3].data[0, 0][data_mask].sum()) * header_stokes["PHOTFLAM"], + Stokes_hdul["I_STOKES"].data[data_mask].sum() * header_stokes["PHOTFLAM"], + np.sqrt(Stokes_hdul["IQU_COV_MATRIX"].data[0, 0][data_mask].sum()) * header_stokes["PHOTFLAM"], 2, out=int, ), diff --git a/package/lib/fits.py b/package/lib/fits.py index 89177c1..6a8c17a 100755 --- a/package/lib/fits.py +++ b/package/lib/fits.py @@ -101,7 +101,24 @@ def get_obs_data(infiles, data_folder="", compute_flux=False): def save_Stokes( - I_stokes, Q_stokes, U_stokes, Stokes_cov, P, debiased_P, s_P, s_P_P, PA, s_PA, s_PA_P, header_stokes, data_mask, filename, data_folder="", return_hdul=False + I_stokes, + Q_stokes, + U_stokes, + Stokes_cov, + P, + debiased_P, + s_P, + s_P_P, + PA, + s_PA, + s_PA_P, + header_stokes, + data_mask, + filename, + data_folder="", + return_hdul=False, + flux_data=None, + flux_head=None, ): """ Save computed polarimetry parameters to a single fits file, @@ -194,11 +211,23 @@ def save_Stokes( hdul = fits.HDUList([]) # Add I_stokes as PrimaryHDU - header["datatype"] = ("I_stokes", "type of data stored in the HDU") - I_stokes[(1 - data_mask).astype(bool)] = 0.0 - primary_hdu = fits.PrimaryHDU(data=I_stokes, header=header) - primary_hdu.name = "I_stokes" - hdul.append(primary_hdu) + if flux_data is None: + header["datatype"] = ("I_stokes", "type of data stored in the HDU") + I_stokes[(1 - data_mask).astype(bool)] = 0.0 + primary_hdu = fits.PrimaryHDU(data=I_stokes, header=header) + primary_hdu.name = "I_stokes" + hdul.append(primary_hdu) + else: + flux_head["TELESCOP"], flux_head["INSTRUME"] = header["TELESCOP"], header["INSTRUME"] + header["datatype"] = ("Flux map", "type of data stored in the HDU") + primary_hdu = fits.PrimaryHDU(data=flux_data, header=flux_head) + primary_hdu.name = "Flux map" + hdul.append(primary_hdu) + header["datatype"] = ("I_stokes", "type of data stored in the HDU") + I_stokes[(1 - data_mask).astype(bool)] = 0.0 + image_hdu = fits.ImageHDU(data=I_stokes, header=header) + image_hdu.name = "I_stokes" + hdul.append(image_hdu) # Add Q, U, Stokes_cov, P, s_P, PA, s_PA to the HDUList for data, name in [ From dea36b25358c82fa42b5ef872f482ba67040d501 Mon Sep 17 00:00:00 2001 From: Thibault Barnouin Date: Wed, 19 Mar 2025 16:38:03 +0100 Subject: [PATCH 02/33] change histogram binning to numpy function --- package/lib/background.py | 32 +++++++++----------------------- 1 file changed, 9 insertions(+), 23 deletions(-) diff --git a/package/lib/background.py b/package/lib/background.py index 99a09f2..96509d2 100755 --- a/package/lib/background.py +++ b/package/lib/background.py @@ -278,7 +278,7 @@ def bkg_fit(data, error, mask, headers, subtract_error=True, display=False, save return n_data_array, n_error_array, headers, background -def bkg_hist(data, error, mask, headers, sub_type=None, subtract_error=True, display=False, savename=None, plots_folder=""): +def bkg_hist(data, error, mask, headers, n_bins=None, subtract_error=True, display=False, savename=None, plots_folder=""): """ ---------- Inputs: @@ -333,29 +333,15 @@ def bkg_hist(data, error, mask, headers, sub_type=None, subtract_error=True, dis for i, image in enumerate(data): # Compute the Count-rate histogram for the image n_mask = np.logical_and(mask, image > 0.0) - if sub_type is not None: - if isinstance(sub_type, int): - n_bins = sub_type - elif sub_type.lower() in ["square-root", "squareroot", "sqrt"]: - n_bins = np.fix(np.sqrt(image[n_mask].size)).astype(int) # Square-root - elif sub_type.lower() in ["sturges"]: - n_bins = np.ceil(np.log2(image[n_mask].size)).astype(int) + 1 # Sturges - elif sub_type.lower() in ["rice"]: - n_bins = 2 * np.fix(np.power(image[n_mask].size, 1 / 3)).astype(int) # Rice - elif sub_type.lower() in ["freedman-diaconis", "freedmandiaconis", "freedman", "diaconis"]: - n_bins = np.fix( - (image[n_mask].max() - image[n_mask].min()) - / (2 * np.subtract(*np.percentile(image[n_mask], [75, 25])) / np.power(image[n_mask].size, 1 / 3)) - ).astype(int) # Freedman-Diaconis - else: # Fallback - n_bins = np.fix((image[n_mask].max() - image[n_mask].min()) / (3.5 * image[n_mask].std() / np.power(image[n_mask].size, 1 / 3))).astype( - int - ) # Scott - else: # Default statistic - n_bins = np.fix((image[n_mask].max() - image[n_mask].min()) / (3.5 * image[n_mask].std() / np.power(image[n_mask].size, 1 / 3))).astype( - int - ) # Scott + if not isinstance(n_bins, int) and n_bins not in ["auto", "fd", "doane", "scott", "stone", "rice", "sturges", "sqrt"]: + match n_bins.lower(): + case "square-root" | "squareroot": + n_bins = "sqrt" + case "freedman-diaconis" | "freedmandiaconis": + n_bins = "fd" + case _: + n_bins = "scott" hist, bin_edges = np.histogram(np.log(image[n_mask]), bins=n_bins) histograms.append(hist) binning.append(np.exp(bin_centers(bin_edges))) From 02f8d5213c84b9f76ea7ffebfb9d1b47b74b8fa9 Mon Sep 17 00:00:00 2001 From: Thibault Barnouin Date: Wed, 19 Mar 2025 16:38:36 +0100 Subject: [PATCH 03/33] save raw flux in fits file and display --- package/FOC_reduction.py | 41 +++--- package/lib/fits.py | 26 ++-- package/lib/plots.py | 279 +++++++++++++++++++++++++-------------- package/lib/reduction.py | 22 ++- 4 files changed, 230 insertions(+), 138 deletions(-) diff --git a/package/FOC_reduction.py b/package/FOC_reduction.py index eb73dc9..899c9e6 100755 --- a/package/FOC_reduction.py +++ b/package/FOC_reduction.py @@ -20,6 +20,7 @@ import lib.reduction as proj_red # Functions used in reduction pipeline import numpy as np from lib.utils import princ_angle, sci_not from matplotlib.colors import LogNorm +from astropy.wcs import WCS def main(target=None, proposal_id=None, infiles=None, output_dir="./data", crop=False, interactive=False): @@ -40,12 +41,12 @@ def main(target=None, proposal_id=None, infiles=None, output_dir="./data", crop= display_crop = False # Background estimation - error_sub_type = "scott" # sqrt, sturges, rice, freedman-diaconis, scott (default) or shape (example (51, 51)) + error_sub_type = "scott" # auto, sqrt, sturges, rice, freedman-diaconis, scott (default) or shape (example (51, 51)) subtract_error = 3.0 - display_bkg = True + display_bkg = False # Data binning - pxsize = 40 + pxsize = 4 pxscale = "px" # pixel, arcsec or full rebin_operation = "sum" # sum or average @@ -63,10 +64,11 @@ def main(target=None, proposal_id=None, infiles=None, output_dir="./data", crop= smoothing_scale = "px" # pixel or arcsec # Rotation + rotate_data = False rotate_North = True # Polarization map output - P_cut = 5 # if >=1.0 cut on the signal-to-noise else cut on the confidence level in Q, U + P_cut = 3 # if >=1.0 cut on the signal-to-noise else cut on the confidence level in Q, U SNRi_cut = 1.0 # I measurments with SNR>30, which implies an uncertainty in P of 4.7%. flux_lim = None # lowest and highest flux displayed on plot, defaults to bkg and maximum in cut if None scale_vec = 3 @@ -142,7 +144,7 @@ def main(target=None, proposal_id=None, infiles=None, output_dir="./data", crop= ) # Rotate data to have same orientation - rotate_data = np.unique([np.round(float(head["ORIENTAT"]), 3) for head in headers]).size != 1 + rotate_data = rotate_data or np.unique([np.round(float(head["ORIENTAT"]), 3) for head in headers]).size != 1 if rotate_data: ang = np.mean([head["ORIENTAT"] for head in headers]) for head in headers: @@ -158,7 +160,6 @@ def main(target=None, proposal_id=None, infiles=None, output_dir="./data", crop= vmin=data_array[data_array > 0.0].min() * headers[0]["photflam"], vmax=data_array[data_array > 0.0].max() * headers[0]["photflam"] ), ) - # Align and rescale images with oversampling. data_array, error_array, headers, data_mask, shifts, error_shifts = proj_red.align_data( data_array, @@ -243,6 +244,8 @@ def main(target=None, proposal_id=None, infiles=None, output_dir="./data", crop= ) flux_data, flux_error, flux_mask, flux_head = proj_red.rotate_data(np.array([flux_data]), np.array([flux_error]), flux_mask, [flux_head]) flux_data, flux_error, flux_head = flux_data[0], flux_error[0], flux_head[0] + elif not rotate_data: + figname += "_noROT" # Compute polarimetric parameters (polarization degree and angle). P, debiased_P, s_P, s_P_P, PA, s_PA, s_PA_P = proj_red.compute_pol(I_stokes, Q_stokes, U_stokes, Stokes_cov, header_stokes) @@ -250,7 +253,7 @@ def main(target=None, proposal_id=None, infiles=None, output_dir="./data", crop= # Step 4: # Save image to FITS. - figname = "_".join([figname, figtype]) if figtype != "" else figname + savename = "_".join([figname, figtype]) if figtype != "" else figname Stokes_hdul = proj_fits.save_Stokes( I_stokes, Q_stokes, @@ -265,18 +268,18 @@ def main(target=None, proposal_id=None, infiles=None, output_dir="./data", crop= s_PA_P, header_stokes, data_mask, - figname, + savename, data_folder=data_folder, return_hdul=True, - flux_data=flux_data, + flux_data=np.array([flux_data, flux_error, flux_mask]), flux_head=flux_head, ) - outfiles.append("/".join([data_folder, Stokes_hdul["I_STOKES"].header["FILENAME"] + ".fits"])) + outfiles.append("/".join([data_folder, Stokes_hdul[0].header["FILENAME"] + ".fits"])) # Step 5: # crop to desired region of interest (roi) if crop: - figname += "_crop" + savename += "_crop" stokescrop = proj_plots.crop_Stokes(deepcopy(Stokes_hdul), norm=LogNorm()) stokescrop.crop() stokescrop.write_to("/".join([data_folder, figname + ".fits"])) @@ -286,10 +289,10 @@ def main(target=None, proposal_id=None, infiles=None, output_dir="./data", crop= data_mask = Stokes_hdul["data_mask"].data.astype(bool) print( "F_int({0:.0f} Angs) = ({1} ± {2})e{3} ergs.cm^-2.s^-1.Angs^-1".format( - header_stokes["PHOTPLAM"], + flux_head["PHOTPLAM"], *sci_not( - Stokes_hdul["I_STOKES"].data[data_mask].sum() * header_stokes["PHOTFLAM"], - np.sqrt(Stokes_hdul["IQU_COV_MATRIX"].data[0, 0][data_mask].sum()) * header_stokes["PHOTFLAM"], + flux_data[flux_mask].sum() * flux_head["PHOTFLAM"], + np.sqrt(np.sum(flux_error[flux_mask] ** 2)) * flux_head["PHOTFLAM"], 2, out=int, ), @@ -315,12 +318,12 @@ def main(target=None, proposal_id=None, infiles=None, output_dir="./data", crop= flux_lim=flux_lim, step_vec=step_vec, scale_vec=scale_vec, - savename="_".join([figname]), + savename=figname, plots_folder=plots_folder, ) for figtype, figsuffix in zip( - ["Intensity", "Pol_flux", "Pol_deg", "Pol_ang", "I_err", "P_err", "SNRi", "SNRp", "confp"], - ["I", "P_flux", "P", "PA", "I_err", "P_err", "SNRi", "SNRp", "confP"], + ["FluxDensity", "Intensity", "Pol_flux", "Pol_deg", "Pol_ang", "I_err", "P_err", "SNRi", "SNRp", "confp"], + ["F", "I", "P_flux", "P", "PA", "I_err", "P_err", "SNRi", "SNRp", "confP"], ): try: proj_plots.polarization_map( @@ -331,7 +334,7 @@ def main(target=None, proposal_id=None, infiles=None, output_dir="./data", crop= flux_lim=flux_lim, step_vec=step_vec, scale_vec=scale_vec, - savename="_".join([figname, figsuffix]), + savename="_".join([savename, figsuffix]), plots_folder=plots_folder, display=figtype, ) @@ -339,7 +342,7 @@ def main(target=None, proposal_id=None, infiles=None, output_dir="./data", crop= pass elif not interactive: proj_plots.polarization_map( - deepcopy(Stokes_hdul), data_mask, P_cut=P_cut, SNRi_cut=SNRi_cut, savename=figname, plots_folder=plots_folder, display="integrate" + deepcopy(Stokes_hdul), data_mask, P_cut=P_cut, SNRi_cut=SNRi_cut, savename=savename, plots_folder=plots_folder, display="integrate" ) elif pxscale.lower() not in ["full", "integrate"]: proj_plots.pol_map(Stokes_hdul, P_cut=P_cut, SNRi_cut=SNRi_cut, step_vec=step_vec, scale_vec=scale_vec, flux_lim=flux_lim) diff --git a/package/lib/fits.py b/package/lib/fits.py index 6a8c17a..feb828a 100755 --- a/package/lib/fits.py +++ b/package/lib/fits.py @@ -47,9 +47,9 @@ def get_obs_data(infiles, data_folder="", compute_flux=False): wcs_array.append(WCS(header=f[0].header, fobj=f).celestial) f.flush() # Save pixel area for flux density computation - if headers[i]["PXFORMT"] == "NORMAL": + if "PXFORMT" in headers[i].keys() and headers[i]["PXFORMT"] == "NORMAL": headers[i]["PXAREA"] = 1.96e-4 # 14x14 milliarcsec squared pixel area in arcsec^2 - elif headers[i]["PXFORMT"] == "ZOOM": + elif "PXFORMT" in headers[i].keys() and headers[i]["PXFORMT"] == "ZOOM": headers[i]["PXAREA"] = 4.06e-4 # 29x14 milliarcsec squared pixel area in arcsec^2 else: headers[i]["PXAREA"] = 1.0 # unknown default to 1 arcsec^2 @@ -85,10 +85,10 @@ def get_obs_data(infiles, data_folder="", compute_flux=False): # force WCS for POL60 to have same pixel size as POL0 and POL120 is_pol60 = np.array([head["filtnam1"].lower() == "pol60" for head in headers], dtype=bool) cdelt = np.round(np.array([WCS(head).wcs.cdelt[:2] for head in headers]), 10) - if np.unique(cdelt[np.logical_not(is_pol60)], axis=0).size != 2: + if np.any(is_pol60) and np.unique(cdelt[np.logical_not(is_pol60)], axis=0).size != 2: print(np.unique(cdelt[np.logical_not(is_pol60)], axis=0)) raise ValueError("Not all images have same pixel size") - else: + elif np.any(is_pol60): for i in np.arange(len(headers))[is_pol60]: headers[i]["cdelt1"], headers[i]["cdelt2"] = np.unique(cdelt[np.logical_not(is_pol60)], axis=0)[0] @@ -174,7 +174,7 @@ def save_Stokes( header["PROPOSID"] = (header_stokes["PROPOSID"], "PEP proposal identifier for observation") header["TARGNAME"] = (header_stokes["TARGNAME"], "Target name") header["ORIENTAT"] = (header_stokes["ORIENTAT"], "Angle between North and the y-axis of the image") - header["FILENAME"] = (filename, "ORIGINAL FILENAME") + header["FILENAME"] = (filename, "Original filename") header["BKG_TYPE"] = (header_stokes["BKG_TYPE"], "Bkg estimation method used during reduction") header["BKG_SUB"] = (header_stokes["BKG_SUB"], "Amount of bkg subtracted from images") header["SMOOTH"] = (header_stokes["SMOOTH"] if "SMOOTH" in list(header_stokes.keys()) else "None", "Smoothing method used during reduction") @@ -210,7 +210,7 @@ def save_Stokes( # Create HDUList object hdul = fits.HDUList([]) - # Add I_stokes as PrimaryHDU + # Add Flux density as PrimaryHDU if flux_data is None: header["datatype"] = ("I_stokes", "type of data stored in the HDU") I_stokes[(1 - data_mask).astype(bool)] = 0.0 @@ -218,10 +218,16 @@ def save_Stokes( primary_hdu.name = "I_stokes" hdul.append(primary_hdu) else: - flux_head["TELESCOP"], flux_head["INSTRUME"] = header["TELESCOP"], header["INSTRUME"] - header["datatype"] = ("Flux map", "type of data stored in the HDU") - primary_hdu = fits.PrimaryHDU(data=flux_data, header=flux_head) - primary_hdu.name = "Flux map" + flux_head["FILENAME"] = header["FILENAME"] + head = WCS(flux_head).deepcopy().to_header() + for key in [key for key in header.keys() if key not in ["SMOOTH", "SAMPLING"]]: + try: + head[key] = flux_head[key] + except KeyError: + head[key] = header[key] + header["datatype"] = ("Flux_density", "type of data stored in the HDU") + primary_hdu = fits.PrimaryHDU(data=flux_data, header=head) + primary_hdu.name = "Flux_density" hdul.append(primary_hdu) header["datatype"] = ("I_stokes", "type of data stored in the HDU") I_stokes[(1 - data_mask).astype(bool)] = 0.0 diff --git a/package/lib/plots.py b/package/lib/plots.py index 40ec567..a3f51ad 100755 --- a/package/lib/plots.py +++ b/package/lib/plots.py @@ -58,7 +58,6 @@ from mpl_toolkits.axes_grid1.anchored_artists import ( AnchoredDirectionArrows, AnchoredSizeBar, ) -from scipy.ndimage import zoom as sc_zoom try: from .utils import PCconf, princ_angle, rot2D, sci_not @@ -166,26 +165,29 @@ def plot_obs( ) ax_curr.plot([x, x], [0, data.shape[0] - 1], "--", lw=2, color="g", alpha=0.85) ax_curr.plot([0, data.shape[1] - 1], [y, y], "--", lw=2, color="g", alpha=0.85) + # position of centroid + ax_curr.plot( + [data.shape[1] / 2, data.shape[1] / 2], + [0, data.shape[0] - 1], + "--", + lw=2, + color="b", + alpha=0.85, + ) + ax_curr.plot( + [0, data.shape[1] - 1], + [data.shape[0] / 2, data.shape[0] / 2], + "--", + lw=2, + color="b", + alpha=0.85, + ) + cr_x, cr_y = head["CRPIX1"], head["CRPIX2"] + # Plot WCS reference point + ax_curr.plot([cr_x], [cr_y], "+", lw=2, color="r", alpha=0.85) if rectangle is not None: x, y, width, height, angle, color = rectangle[i] ax_curr.add_patch(Rectangle((x, y), width, height, angle=angle, edgecolor=color, fill=False)) - # position of centroid - ax_curr.plot( - [data.shape[1] / 2, data.shape[1] / 2], - [0, data.shape[0] - 1], - "--", - lw=2, - color="b", - alpha=0.85, - ) - ax_curr.plot( - [0, data.shape[1] - 1], - [data.shape[0] / 2, data.shape[0] / 2], - "--", - lw=2, - color="b", - alpha=0.85, - ) ax_curr.annotate( instr + ":" + rootname, color="white", @@ -269,7 +271,7 @@ def plot_Stokes(Stokes, savename=None, plots_folder=""): for dataset in [stkI, stkQ, stkU]: dataset[np.logical_not(data_mask)] = np.nan - wcs = WCS(Stokes[0]).deepcopy() + wcs = WCS(Stokes["I_STOKES"]).deepcopy() # Plot figure plt.rcParams.update({"font.size": 14}) @@ -373,6 +375,9 @@ def polarization_map( The figure and ax created for interactive contour maps. """ # Get data + flux = Stokes[0].data[0].copy() * Stokes[0].header["PHOTFLAM"] + flux_error = Stokes[0].data[1].copy() * Stokes[0].header["PHOTFLAM"] + flux_mask = Stokes[0].data[2].astype(bool).copy() stkI = Stokes["I_stokes"].data.copy() stkQ = Stokes["Q_stokes"].data.copy() stkU = Stokes["U_stokes"].data.copy() @@ -387,6 +392,20 @@ def polarization_map( data_mask = np.zeros(stkI.shape).astype(bool) data_mask[stkI > 0.0] = True + wcs = WCS(Stokes["I_STOKES"]).deepcopy() + pivot_wav = Stokes["I_STOKES"].header["photplam"] + convert_flux = Stokes["I_STOKES"].header["photflam"] + + # Get integrated flux values from sum + I_diluted = stkI[data_mask].sum() * convert_flux + I_diluted_err = np.sqrt(np.sum(stk_cov[0, 0][data_mask])) * convert_flux + + # Get integrated polarization values from header + P_diluted = Stokes["I_STOKES"].header["P_int"] + P_diluted_err = Stokes["I_STOKES"].header["sP_int"] + PA_diluted = Stokes["I_STOKES"].header["PA_int"] + PA_diluted_err = Stokes["I_STOKES"].header["sPA_int"] + # Compute confidence level map QN, UN, QN_ERR, UN_ERR = np.full((4, stkI.shape[0], stkI.shape[1]), np.nan) for nflux, sflux in zip( @@ -402,12 +421,8 @@ def polarization_map( for j in range(3): stk_cov[i][j][np.logical_not(data_mask)] = np.nan - wcs = WCS(Stokes[0]).deepcopy() - pivot_wav = Stokes[0].header["photplam"] - convert_flux = Stokes[0].header["photflam"] - # Plot Stokes parameters map - if display is None or display.lower() in ["default"]: + if display is None or display.lower() in ["pol", "polarization", "polarisation", "pol_deg", "p"]: plot_Stokes(Stokes, savename=savename, plots_folder=plots_folder) # Compute SNR and apply cuts @@ -451,7 +466,7 @@ def polarization_map( fig = plt.figure(figsize=(7 * ratiox, 7 * ratioy), layout="constrained") if ax is None: ax = fig.add_subplot(111, projection=wcs) - ax.set(aspect="equal", fc="k") + ax.set(aspect="equal", fc="k", xlim=(0, stkI.shape[1]), ylim=(0, stkI.shape[0])) # fig.subplots_adjust(hspace=0, wspace=0, left=0.102, right=1.02) # ax.coords.grid(True, color='white', ls='dotted', alpha=0.5) @@ -501,7 +516,30 @@ def polarization_map( ax.set_facecolor("white") font_color = "black" - if display.lower() in ["i", "intensity"]: + if display.lower() in ["f", "flux", "fluxdensity"]: + # If no display selected, show intensity map + display = "f" + if flux_lim is not None: + vmin, vmax = flux_lim + else: + vmin, vmax = np.max(flux[flux > 0.0]) / 2e3, np.max(flux[flux > 0.0]) + imflux, cr = flux.copy(), WCS(Stokes[0].header).wcs.crpix.astype(int) + imflux[cr[1] - 1 : cr[1] + 2, cr[0] - 1 : cr[0] + 2] = np.nan + im = ax.imshow( + imflux, + transform=ax.get_transform(WCS(Stokes[0].header).celestial), + norm=LogNorm(vmin, vmax), + aspect="equal", + cmap=kwargs["cmap"], + alpha=1.0, + ) + fig.colorbar(im, ax=ax, aspect=50, shrink=0.60, pad=0.025, label=r"$F_{\lambda}$ [$ergs \cdot cm^{-2} \cdot s^{-1} \cdot \AA$]") + levelsF = np.array([0.8, 2.0, 5.0, 10.0, 20.0, 50.0]) / 100.0 * vmax + print("Flux density contour levels : ", levelsF) + ax.contour(flux, levels=levelsF, transform=ax.get_transform(WCS(Stokes[0].header).celestial), colors="grey", linewidths=0.5) + ax.plot(*WCS(Stokes[1]).wcs.crpix, "g+") + I_diluted, I_diluted_err = np.sum(flux[flux_mask]), np.sqrt(np.sum(flux_error[flux_mask] ** 2)) + elif display.lower() in ["i", "intensity"]: # If no display selected, show intensity map display = "i" if flux_lim is not None: @@ -512,9 +550,12 @@ def polarization_map( vmin, vmax = 1.0 / 2.0 * np.median(np.sqrt(stk_cov[0, 0][stkI > 0.0]) * convert_flux), np.max(stkI[stkI > 0.0] * convert_flux) im = ax.imshow(stkI * convert_flux, norm=LogNorm(vmin, vmax), aspect="equal", cmap=kwargs["cmap"], alpha=1.0) fig.colorbar(im, ax=ax, aspect=50, shrink=0.60, pad=0.025, label=r"$F_{\lambda}$ [$ergs \cdot cm^{-2} \cdot s^{-1} \cdot \AA^{-1}$]") - levelsI = np.array([0.8, 2.0, 5.0, 10.0, 20.0, 50.0]) / 100.0 * vmax - print("Flux density contour levels : ", levelsI) - ax.contour(stkI * convert_flux, levels=levelsI, colors="grey", linewidths=0.5) + # levelsI = np.array([0.8, 2.0, 5.0, 10.0, 20.0, 50.0]) / 100.0 * vmax + # print("Stokes I contour levels : ", levelsI) + # ax.contour(stkI * convert_flux, levels=levelsI, colors="grey", linewidths=0.5) + levelsF = np.array([0.8, 2.0, 5.0, 10.0, 20.0, 50.0]) / 100.0 * np.max(flux[flux > 0.0]) + print("Flux density contour levels : ", levelsF) + ax.contour(flux, levels=levelsF, transform=ax.get_transform(WCS(Stokes[0].header).celestial), colors="grey", linewidths=0.5) elif display.lower() in ["pf", "pol_flux"]: # Display polarization flux display = "pf" @@ -604,22 +645,18 @@ def polarization_map( # Defaults to intensity map if flux_lim is not None: vmin, vmax = flux_lim - elif mask.sum() > 0.0: - vmin, vmax = 1.0 * np.mean(np.sqrt(stk_cov[0, 0][mask]) * convert_flux), np.max(stkI[stkI > 0.0] * convert_flux) else: - vmin, vmax = 1.0 * np.mean(np.sqrt(stk_cov[0, 0][stkI > 0.0]) * convert_flux), np.max(stkI[stkI > 0.0] * convert_flux) - im = ax.imshow(stkI * convert_flux, norm=LogNorm(vmin, vmax), aspect="equal", cmap=kwargs["cmap"], alpha=1.0) + vmin, vmax = np.max(flux[flux > 0.0] * convert_flux) / 2e3, np.max(flux[flux > 0.0] * convert_flux) + im = ax.imshow( + flux * Stokes[0].header["PHOTFLAM"], + transform=ax.get_transform(WCS(Stokes[0].header).celestial), + norm=LogNorm(vmin, vmax), + aspect="equal", + cmap=kwargs["cmap"], + alpha=1.0, + ) fig.colorbar(im, ax=ax, aspect=50, shrink=0.60, pad=0.025, label=r"$F_{\lambda}$ [$ergs \cdot cm^{-2} \cdot s^{-1} \cdot \AA$]") - - # Get integrated flux values from sum - I_diluted = stkI[data_mask].sum() - I_diluted_err = np.sqrt(np.sum(stk_cov[0, 0][data_mask])) - - # Get integrated polarization values from header - P_diluted = Stokes[0].header["P_int"] - P_diluted_err = Stokes[0].header["sP_int"] - PA_diluted = Stokes[0].header["PA_int"] - PA_diluted_err = Stokes[0].header["sPA_int"] + I_diluted, I_diluted_err = np.sum(flux[flux_mask]), np.sqrt(np.sum(flux_error[flux_mask] ** 2)) plt.rcParams.update({"font.size": 11}) px_size = wcs.wcs.get_cdelt()[0] * 3600.0 @@ -648,12 +685,12 @@ def polarization_map( back_length=0.0, head_length=7.0, head_width=7.0, - angle=-Stokes[0].header["orientat"], + angle=-Stokes["I_STOKES"].header["orientat"], text_props={"ec": "k", "fc": font_color, "alpha": 1, "lw": 0.5}, arrow_props={"ec": "k", "fc": font_color, "alpha": 1, "lw": 1}, ) - if display.lower() in ["i", "s_i", "snri", "pf", "p", "pa", "s_p", "snrp", "confp"] and step_vec != 0: + if display.lower() in ["f", "i", "s_i", "snri", "pf", "p", "pa", "s_p", "snrp", "confp"] and step_vec != 0: if scale_vec == -1: poldata[np.isfinite(poldata)] = 1.0 / 2.0 step_vec = 1 @@ -701,7 +738,7 @@ def polarization_map( ax.annotate( r"$F_{{\lambda}}^{{int}}$({0:.0f} $\AA$) = {1} $ergs \cdot cm^{{-2}} \cdot s^{{-1}} \cdot \AA^{{-1}}$".format( pivot_wav, - sci_not(I_diluted * convert_flux, I_diluted_err * convert_flux, 2), + sci_not(I_diluted, I_diluted_err, 2), ) + "\n" + r"$P^{{int}}$ = {0:.1f} $\pm$ {1:.1f} %".format(P_diluted * 100.0, P_diluted_err * 100.0) @@ -721,7 +758,7 @@ def polarization_map( ax.annotate( r"$F_{{\lambda}}^{{int}}$({0:.0f} $\AA$) = {1} $ergs \cdot cm^{{-2}} \cdot s^{{-1}} \cdot \AA^{{-1}}$".format( pivot_wav, - sci_not(I_diluted * convert_flux, I_diluted_err * convert_flux, 2), + sci_not(I_diluted, I_diluted_err, 2), ), color="white", xy=(0.01, 1.00), @@ -1373,6 +1410,8 @@ class overplot_chandra(align_maps): other_data = deepcopy(self.other_data) other_wcs = self.other_wcs.deepcopy() if zoom != 1: + from scipy.ndimage import zoom as sc_zoom + other_data = sc_zoom(other_data, zoom) other_wcs.wcs.crpix *= zoom other_wcs.wcs.cdelt /= zoom @@ -2947,9 +2986,9 @@ class pol_map(object): self.conf = PCconf(self.QN, self.UN, self.QN_ERR, self.UN_ERR) # Get data - self.targ = self.Stokes[0].header["targname"] - self.pivot_wav = self.Stokes[0].header["photplam"] - self.map_convert = self.Stokes[0].header["photflam"] + self.targ = self.Stokes["I_STOKES"].header["targname"] + self.pivot_wav = self.Stokes["I_STOKES"].header["photplam"] + self.map_convert = self.Stokes["I_STOKES"].header["photflam"] # Create figure plt.rcParams.update({"font.size": 10}) @@ -3062,7 +3101,7 @@ class pol_map(object): def select_roi(event): if self.data is None: - self.data = self.Stokes[0].data + self.data = self.Stokes["I_STOKES"].data if self.selected: self.selected = False self.region = deepcopy(self.select_instance.mask.astype(bool)) @@ -3111,7 +3150,7 @@ class pol_map(object): def select_aperture(event): if self.data is None: - self.data = self.Stokes[0].data + self.data = self.Stokes["I_STOKES"].data if self.selected: self.selected = False self.select_instance.update_mask() @@ -3200,7 +3239,7 @@ class pol_map(object): def select_slit(event): if self.data is None: - self.data = self.Stokes[0].data + self.data = self.Stokes["I_STOKES"].data if self.selected: self.selected = False self.select_instance.update_mask() @@ -3534,7 +3573,19 @@ class pol_map(object): @property def wcs(self): - return WCS(self.Stokes[0].header).celestial.deepcopy() + return WCS(self.Stokes["I_STOKES"].header).celestial.deepcopy() + + @property + def Flux(self): + return self.Stokes[0].data[0] * self.Stokes[0].header["PHOTFLAM"] + + @property + def Flux_err(self): + return self.Stokes[0].data[1] * self.Stokes[0].header["PHOTFLAM"] + + @property + def Flux_mask(self): + return self.Stokes[0].data[2].astype(bool) @property def I(self): @@ -3598,7 +3649,7 @@ class pol_map(object): @property def data_mask(self): - return self.Stokes["DATA_MASK"].data + return self.Stokes["DATA_MASK"].data.astype(bool) def set_data_mask(self, mask): self.Stokes[np.argmax([self.Stokes[i].header["datatype"] == "Data_mask" for i in range(len(self.Stokes))])].data = mask.astype(float) @@ -3682,7 +3733,7 @@ class pol_map(object): back_length=0.0, head_length=10.0, head_width=10.0, - angle=-self.Stokes[0].header["orientat"], + angle=-self.Stokes["I_STOKES"].header["orientat"], color="white", text_props={"ec": None, "fc": "w", "alpha": 1, "lw": 0.4}, arrow_props={"ec": None, "fc": "w", "alpha": 1, "lw": 1}, @@ -3690,21 +3741,23 @@ class pol_map(object): ax.add_artist(self.north_dir) def display(self, fig=None, ax=None, flux_lim=None): - norm = None - if self.display_selection is None: - self.display_selection = "total_flux" + kwargs = dict([]) if flux_lim is None: flux_lim = self.flux_lim - if self.display_selection.lower() in ["total_flux"]: - self.data = self.I * self.map_convert + if self.display_selection is None or self.display_selection.lower() in ["total_flux"]: + self.data = self.Flux # self.I * self.map_convert if flux_lim is None: vmin, vmax = ( - 1.0 / 2.0 * np.median(self.data[self.data > 0.0]), - np.max(self.data[self.data > 0.0]), + 1.0 / 2.0 * np.median(self.I[self.I > 0.0] * self.map_convert), + np.max(self.I[self.I > 0.0] * self.map_convert), ) else: vmin, vmax = flux_lim - norm = LogNorm(vmin, vmax) + kwargs["norm"] = LogNorm(vmin, vmax) + if ax is None: + kwargs["transform"] = self.ax.get_transform(WCS(self.Stokes[0].header).celestial) + else: + kwargs["transform"] = ax.get_transform(WCS(self.Stokes[0].header).celestial) label = r"$F_{\lambda}$ [$ergs \cdot cm^{-2} \cdot s^{-1} \cdot \AA^{-1}$]" elif self.display_selection.lower() in ["pol_flux"]: self.data = self.I * self.map_convert * self.P @@ -3715,28 +3768,28 @@ class pol_map(object): ) else: vmin, vmax = flux_lim - norm = LogNorm(vmin, vmax) + kwargs["norm"] = LogNorm(vmin, vmax) label = r"$P \cdot F_{\lambda}$ [$ergs \cdot cm^{-2} \cdot s^{-1} \cdot \AA^{-1}$]" elif self.display_selection.lower() in ["pol_deg"]: self.data = self.P * 100.0 - vmin, vmax = 0.0, np.max(self.data[self.P > self.s_P]) + kwargs["vmin"], kwargs["vmax"] = 0.0, np.max(self.data[self.P > self.s_P]) label = r"$P$ [%]" elif self.display_selection.lower() in ["pol_ang"]: self.data = princ_angle(self.PA) - vmin, vmax = 0, 180.0 + kwargs["vmin"], kwargs["vmax"] = 0, 180.0 label = r"$\theta_{P}$ [°]" elif self.display_selection.lower() in ["snri"]: s_I = np.sqrt(self.IQU_cov[0, 0]) SNRi = np.zeros(self.I.shape) SNRi[s_I > 0.0] = self.I[s_I > 0.0] / s_I[s_I > 0.0] self.data = SNRi - vmin, vmax = 0.0, np.max(self.data[self.data > 0.0]) + kwargs["vmin"], kwargs["vmax"] = 0.0, np.max(self.data[self.data > 0.0]) label = r"$I_{Stokes}/\sigma_{I}$" elif self.display_selection.lower() in ["snrp"]: SNRp = np.zeros(self.P.shape) SNRp[self.s_P > 0.0] = self.P[self.s_P > 0.0] / self.s_P[self.s_P > 0.0] self.data = SNRp - vmin, vmax = 0.0, np.max(self.data[self.data > 0.0]) + kwargs["vmin"], kwargs["vmax"] = 0.0, np.max(self.data[self.data > 0.0]) label = r"$P/\sigma_{P}$" if fig is None: @@ -3747,22 +3800,17 @@ class pol_map(object): self.cbar.remove() if hasattr(self, "im"): self.im.remove() - if norm is not None: - self.im = ax.imshow(self.data, norm=norm, aspect="equal", cmap="inferno") - else: - self.im = ax.imshow(self.data, vmin=vmin, vmax=vmax, aspect="equal", cmap="inferno") + self.im = ax.imshow(self.data, aspect="equal", cmap="inferno", **kwargs) + ax.set(xlim=(0, self.I.shape[1]), ylim=(0, self.I.shape[0])) plt.rcParams.update({"font.size": 14}) self.cbar = fig.colorbar(self.im, ax=ax, aspect=50, shrink=0.75, pad=0.025, label=label) plt.rcParams.update({"font.size": 10}) fig.canvas.draw_idle() return self.im else: - if norm is not None: - im = ax.imshow(self.data, norm=norm, aspect="equal", cmap="inferno") - else: - im = ax.imshow(self.data, vmin=vmin, vmax=vmax, aspect="equal", cmap="inferno") - ax.set_xlim(0, self.data.shape[1]) - ax.set_ylim(0, self.data.shape[0]) + im = ax.imshow(self.data, aspect="equal", cmap="inferno", **kwargs) + ax.set_xlim(0, self.I.shape[1]) + ax.set_ylim(0, self.I.shape[0]) plt.rcParams.update({"font.size": 14}) fig.colorbar(im, ax=ax, aspect=50, shrink=0.75, pad=0.025, label=label) plt.rcParams.update({"font.size": 10}) @@ -3789,24 +3837,24 @@ class pol_map(object): ax = self.ax if hasattr(self, "quiver"): self.quiver.remove() - self.quiver = ax.quiver( - X[:: self.step_vec, :: self.step_vec], - Y[:: self.step_vec, :: self.step_vec], - XY_U[:: self.step_vec, :: self.step_vec], - XY_V[:: self.step_vec, :: self.step_vec], - units="xy", - scale=1.0 / scale_vec, - scale_units="xy", - pivot="mid", - headwidth=0.0, - headlength=0.0, - headaxislength=0.0, - width=0.3, - linewidth=0.6, - color="white", - edgecolor="black", - ) if self.pa_err: + self.quiver = ax.quiver( + X[:: self.step_vec, :: self.step_vec], + Y[:: self.step_vec, :: self.step_vec], + XY_U[:: self.step_vec, :: self.step_vec], + XY_V[:: self.step_vec, :: self.step_vec], + units="xy", + scale=1.0 / scale_vec, + scale_units="xy", + pivot="mid", + headwidth=0.0, + headlength=0.0, + headaxislength=0.0, + width=0.1, + # linewidth=0.6, + color="black", + edgecolor="black", + ) XY_U_err1, XY_V_err1 = ( P_cut * np.cos(np.pi / 2.0 + (self.PA + 3.0 * self.s_PA) * np.pi / 180.0), P_cut * np.sin(np.pi / 2.0 + (self.PA + 3.0 * self.s_PA) * np.pi / 180.0), @@ -3855,6 +3903,25 @@ class pol_map(object): edgecolor="black", ls="dashed", ) + else: + self.quiver = ax.quiver( + X[:: self.step_vec, :: self.step_vec], + Y[:: self.step_vec, :: self.step_vec], + XY_U[:: self.step_vec, :: self.step_vec], + XY_V[:: self.step_vec, :: self.step_vec], + units="xy", + scale=1.0 / scale_vec, + scale_units="xy", + pivot="mid", + headwidth=0.0, + headlength=0.0, + headaxislength=0.0, + width=0.3, + linewidth=0.6, + color="white", + edgecolor="black", + ) + fig.canvas.draw_idle() return self.quiver else: @@ -3926,12 +3993,20 @@ class pol_map(object): str_conf = "" if self.region is None: s_I = np.sqrt(self.IQU_cov[0, 0]) - I_reg = self.I.sum() - I_reg_err = np.sqrt(np.sum(s_I**2)) - P_reg = self.Stokes[0].header["P_int"] - P_reg_err = self.Stokes[0].header["sP_int"] - PA_reg = self.Stokes[0].header["PA_int"] - PA_reg_err = self.Stokes[0].header["sPA_int"] + I_reg = ( + np.sum(self.Flux[self.Flux_mask]) / self.map_convert + if self.display_selection is None or self.display_selection.lower() in ["total_flux"] + else np.sum(self.I[self.data_mask]) + ) + I_reg_err = ( + np.sqrt(np.sum(self.Flux_err[self.Flux_mask] ** 2)) / self.map_convert + if self.display_selection is None or self.display_selection.lower() in ["total_flux"] + else np.sqrt(np.sum(s_I[self.data_mask] ** 2)) + ) + P_reg = self.Stokes["I_STOKES"].header["P_int"] + P_reg_err = self.Stokes["I_STOKES"].header["sP_int"] + PA_reg = self.Stokes["I_STOKES"].header["PA_int"] + PA_reg_err = self.Stokes["I_STOKES"].header["sPA_int"] s_I = np.sqrt(self.IQU_cov[0, 0]) s_Q = np.sqrt(self.IQU_cov[1, 1]) diff --git a/package/lib/reduction.py b/package/lib/reduction.py index 942a27f..a38e7db 100755 --- a/package/lib/reduction.py +++ b/package/lib/reduction.py @@ -52,7 +52,6 @@ from scipy.ndimage import rotate as sc_rotate from scipy.ndimage import shift as sc_shift from scipy.signal import fftconvolve -from .background import bkg_fit, bkg_hist, bkg_mini from .convex_hull import image_hull from .cross_correlation import phase_cross_correlation from .deconvolve import deconvolve_im, gaussian2d, gaussian_psf, zeropad @@ -310,6 +309,8 @@ def crop_array(data_array, headers, error_array=None, data_mask=None, step=5, nu # Update CRPIX value in the associated header curr_wcs = WCS(crop_headers[i]).celestial.deepcopy() curr_wcs.wcs.crpix[:2] = curr_wcs.wcs.crpix[:2] - np.array([v_array[2], v_array[0]]) + curr_wcs.array_shape = crop_array[i].shape + curr_wcs.wcs.set() crop_headers[i].update(curr_wcs.to_header()) crop_headers[i]["naxis1"], crop_headers[i]["naxis2"] = crop_array[i].shape @@ -528,20 +529,22 @@ def get_error( # estimated to less than 3% err_flat = data * 0.03 + from .background import bkg_fit, bkg_hist, bkg_mini + if sub_type is None: n_data_array, c_error_bkg, headers, background = bkg_hist( data, error, mask, headers, subtract_error=subtract_error, display=display, savename=savename, plots_folder=plots_folder ) sub_type, subtract_error = "histogram ", str(int(subtract_error > 0.0)) elif isinstance(sub_type, str): - if sub_type.lower() in ["auto"]: + if sub_type.lower() in ["fit"]: n_data_array, c_error_bkg, headers, background = bkg_fit( data, error, mask, headers, subtract_error=subtract_error, display=display, savename=savename, plots_folder=plots_folder ) sub_type, subtract_error = "histogram fit ", "mean+%.1fsigma" % subtract_error if subtract_error != 0.0 else 0.0 else: n_data_array, c_error_bkg, headers, background = bkg_hist( - data, error, mask, headers, sub_type=sub_type, subtract_error=subtract_error, display=display, savename=savename, plots_folder=plots_folder + data, error, mask, headers, n_bins=sub_type, subtract_error=subtract_error, display=display, savename=savename, plots_folder=plots_folder ) sub_type, subtract_error = "histogram ", "mean+%.1fsigma" % subtract_error if subtract_error != 0.0 else 0.0 elif isinstance(sub_type, tuple): @@ -675,7 +678,9 @@ def rebin_array(data_array, error_array, headers, pxsize=2, scale="px", operatio nw.wcs.cdelt *= Dxy nw.wcs.crpix /= Dxy nw.array_shape = new_shape + nw.wcs.set() new_header["NAXIS1"], new_header["NAXIS2"] = nw.array_shape + new_header["PXAREA"] *= Dxy[0] * Dxy[1] for key, val in nw.to_header().items(): new_header.set(key, val) new_header["SAMPLING"] = (str(pxsize) + scale, "Resampling performed during reduction") @@ -853,7 +858,10 @@ def align_data( new_crpix = np.array([wcs.wcs.crpix for wcs in headers_wcs]) + shifts[:, ::-1] + res_shift[::-1] for i in range(len(headers_wcs)): headers_wcs[i].wcs.crpix = new_crpix[0] + headers_wcs[i].array_shape = (res_shape, res_shape) + headers_wcs[i].wcs.set() headers[i].update(headers_wcs[i].to_header()) + headers[i]["NAXIS1"], headers[i]["NAXIS2"] = res_shape, res_shape data_mask = rescaled_mask.all(axis=0) data_array, error_array, data_mask, headers = crop_array(rescaled_image, headers, rescaled_error, data_mask, null_val=0.01 * background) @@ -1722,13 +1730,11 @@ def rotate_Stokes(I_stokes, Q_stokes, U_stokes, Stokes_cov, data_mask, header_st new_wcs.wcs.pc = np.dot(mrot, new_wcs.wcs.pc) new_wcs.wcs.crpix = np.dot(mrot, new_wcs.wcs.crpix - old_center[::-1]) + new_center[::-1] + new_wcs.array_shape = shape new_wcs.wcs.set() for key, val in new_wcs.to_header().items(): new_header_stokes.set(key, val) - if new_wcs.wcs.pc[0, 0] == 1.0: - new_header_stokes.set("PC1_1", 1.0) - if new_wcs.wcs.pc[1, 1] == 1.0: - new_header_stokes.set("PC2_2", 1.0) + new_header_stokes["NAXIS1"], new_header_stokes["NAXIS2"] = new_wcs.array_shape new_header_stokes["ORIENTAT"] += ang # Nan handling : @@ -1829,9 +1835,11 @@ def rotate_data(data_array, error_array, data_mask, headers): new_wcs = WCS(header).celestial.deepcopy() new_wcs.wcs.pc[:2, :2] = np.dot(mrot, new_wcs.wcs.pc[:2, :2]) new_wcs.wcs.crpix[:2] = np.dot(mrot, new_wcs.wcs.crpix[:2] - old_center[::-1]) + new_center[::-1] + new_wcs.array_shape = shape new_wcs.wcs.set() for key, val in new_wcs.to_header().items(): new_header[key] = val + new_header["NAXIS1"], new_header["NAXIS2"] = new_wcs.array_shape new_header["ORIENTAT"] = np.arccos(new_wcs.celestial.wcs.pc[0, 0]) * 180.0 / np.pi new_header["ROTATE"] = ang new_headers.append(new_header) From 71e9d89091c444d6f28c3a827174ce63d83d7d98 Mon Sep 17 00:00:00 2001 From: Thibault Barnouin Date: Tue, 1 Apr 2025 16:20:40 +0200 Subject: [PATCH 04/33] looking for displacement of WCS in pipeline --- package/FOC_reduction.py | 1 - package/lib/plots.py | 2 +- package/lib/reduction.py | 4 +++- 3 files changed, 4 insertions(+), 3 deletions(-) diff --git a/package/FOC_reduction.py b/package/FOC_reduction.py index 899c9e6..b42f57f 100755 --- a/package/FOC_reduction.py +++ b/package/FOC_reduction.py @@ -20,7 +20,6 @@ import lib.reduction as proj_red # Functions used in reduction pipeline import numpy as np from lib.utils import princ_angle, sci_not from matplotlib.colors import LogNorm -from astropy.wcs import WCS def main(target=None, proposal_id=None, infiles=None, output_dir="./data", crop=False, interactive=False): diff --git a/package/lib/plots.py b/package/lib/plots.py index a3f51ad..a325b51 100755 --- a/package/lib/plots.py +++ b/package/lib/plots.py @@ -549,12 +549,12 @@ def polarization_map( else: vmin, vmax = 1.0 / 2.0 * np.median(np.sqrt(stk_cov[0, 0][stkI > 0.0]) * convert_flux), np.max(stkI[stkI > 0.0] * convert_flux) im = ax.imshow(stkI * convert_flux, norm=LogNorm(vmin, vmax), aspect="equal", cmap=kwargs["cmap"], alpha=1.0) + ax.plot(*WCS(Stokes[1]).wcs.crpix, "g+") fig.colorbar(im, ax=ax, aspect=50, shrink=0.60, pad=0.025, label=r"$F_{\lambda}$ [$ergs \cdot cm^{-2} \cdot s^{-1} \cdot \AA^{-1}$]") # levelsI = np.array([0.8, 2.0, 5.0, 10.0, 20.0, 50.0]) / 100.0 * vmax # print("Stokes I contour levels : ", levelsI) # ax.contour(stkI * convert_flux, levels=levelsI, colors="grey", linewidths=0.5) levelsF = np.array([0.8, 2.0, 5.0, 10.0, 20.0, 50.0]) / 100.0 * np.max(flux[flux > 0.0]) - print("Flux density contour levels : ", levelsF) ax.contour(flux, levels=levelsF, transform=ax.get_transform(WCS(Stokes[0].header).celestial), colors="grey", linewidths=0.5) elif display.lower() in ["pf", "pol_flux"]: # Display polarization flux diff --git a/package/lib/reduction.py b/package/lib/reduction.py index a38e7db..65639fe 100755 --- a/package/lib/reduction.py +++ b/package/lib/reduction.py @@ -647,7 +647,8 @@ def rebin_array(data_array, error_array, headers, pxsize=2, scale="px", operatio pxsize, scale = "", "full" else: raise ValueError("'{0:s}' invalid scale for binning.".format(scale)) - new_shape = np.ceil(min(image.shape / Dxy_arr, key=lambda x: x[0] + x[1])).astype(int) + new_shape_float = min(image.shape / Dxy_arr, key=lambda x: x[0] + x[1]) + new_shape = np.ceil(new_shape_float).astype(int) for i, (image, error, header) in enumerate(list(zip(data_array, error_array, headers))): # Get current pixel size @@ -676,6 +677,7 @@ def rebin_array(data_array, error_array, headers, pxsize=2, scale="px", operatio # Update header nw = w.deepcopy() nw.wcs.cdelt *= Dxy + # nw.wcs.crpix += np.abs(new_shape_float - new_shape) * np.array(new_shape) / Dxy nw.wcs.crpix /= Dxy nw.array_shape = new_shape nw.wcs.set() From b71aaa37e657011d0c23cccf63a4b2bad1dc8d26 Mon Sep 17 00:00:00 2001 From: Thibault Barnouin Date: Tue, 11 Mar 2025 16:14:14 +0100 Subject: [PATCH 05/33] small improvments to ConfCenter et emission_center --- package/lib/utils.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/package/lib/utils.py b/package/lib/utils.py index ee9fc5b..b76b063 100755 --- a/package/lib/utils.py +++ b/package/lib/utils.py @@ -109,7 +109,7 @@ def CenterConf(mask, PA, sPA): conf[np.isfinite(PA)] = 0.5 * np.exp(-0.5 * chi2[np.isfinite(PA)]) c0 = np.unravel_index(np.argmax(conf), conf.shape)[::-1] - result = minimize(chisq, c0, bounds=[(0, PA.shape[1]), (0.0, PA.shape[0])], method="trust-constr") + result = minimize(chisq, c0, bounds=[(0, PA.shape[1]), (0.0, PA.shape[0])], method="TNC") if result.success: print("Center of emission found: reduced chi_squared {0:.2f}/{1:d}={2:.2f}".format(chisq(result.x) * (Nobs - 2), Nobs - 2, chisq(result.x))) else: From 84bbdb9634d8ddf8f114c4c377eaebba37a5b5f9 Mon Sep 17 00:00:00 2001 From: Thibault Barnouin Date: Fri, 14 Mar 2025 14:20:38 +0100 Subject: [PATCH 06/33] rollback CenterConf to simple CDF of chi2 for 2 dof --- package/lib/utils.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/package/lib/utils.py b/package/lib/utils.py index b76b063..ee9fc5b 100755 --- a/package/lib/utils.py +++ b/package/lib/utils.py @@ -109,7 +109,7 @@ def CenterConf(mask, PA, sPA): conf[np.isfinite(PA)] = 0.5 * np.exp(-0.5 * chi2[np.isfinite(PA)]) c0 = np.unravel_index(np.argmax(conf), conf.shape)[::-1] - result = minimize(chisq, c0, bounds=[(0, PA.shape[1]), (0.0, PA.shape[0])], method="TNC") + result = minimize(chisq, c0, bounds=[(0, PA.shape[1]), (0.0, PA.shape[0])], method="trust-constr") if result.success: print("Center of emission found: reduced chi_squared {0:.2f}/{1:d}={2:.2f}".format(chisq(result.x) * (Nobs - 2), Nobs - 2, chisq(result.x))) else: From 075dc4fd1561ccff1b997ec1385d09038a12d105 Mon Sep 17 00:00:00 2001 From: Thibault Barnouin Date: Tue, 1 Apr 2025 17:02:31 +0200 Subject: [PATCH 07/33] fix WCS computation, cdelt should not be sorted --- package/lib/fits.py | 6 +++--- package/lib/utils.py | 2 +- 2 files changed, 4 insertions(+), 4 deletions(-) diff --git a/package/lib/fits.py b/package/lib/fits.py index feb828a..4ea9a47 100755 --- a/package/lib/fits.py +++ b/package/lib/fits.py @@ -16,7 +16,7 @@ from astropy.io import fits from astropy.wcs import WCS from .convex_hull import clean_ROI -from .utils import wcs_PA +from .utils import wcs_PA, princ_angle def get_obs_data(infiles, data_folder="", compute_flux=False): @@ -72,13 +72,13 @@ def get_obs_data(infiles, data_folder="", compute_flux=False): for key in keys: header.remove(key, ignore_missing=True) new_cdelt = np.linalg.eigvals(wcs.wcs.cd) - new_cdelt.sort() + # new_cdelt.sort() new_wcs.wcs.pc = wcs.wcs.cd.dot(np.diag(1.0 / new_cdelt)) new_wcs.wcs.cdelt = new_cdelt for key, val in new_wcs.to_header().items(): header[key] = val try: - _ = header["ORIENTAT"] + header["ORIENTAT"] = princ_angle(float(header["ORIENTAT"])) except KeyError: header["ORIENTAT"] = wcs_PA(new_wcs.wcs.pc[1, 0], np.diag(new_wcs.wcs.pc).mean()) diff --git a/package/lib/utils.py b/package/lib/utils.py index ee9fc5b..9e36dbc 100755 --- a/package/lib/utils.py +++ b/package/lib/utils.py @@ -177,4 +177,4 @@ def wcs_PA(PC21, PC22): orient = np.arccos(PC22) * 180.0 / np.pi elif (abs(PC21) < abs(PC22)) and (PC22 < 0): orient = -np.arccos(PC22) * 180.0 / np.pi - return orient + return princ_angle(orient) From 6387e2352578e5bcea27e9b73d59b49c97c2f0e3 Mon Sep 17 00:00:00 2001 From: Thibault Barnouin Date: Fri, 14 Mar 2025 14:30:30 +0100 Subject: [PATCH 08/33] Save the raw total flux image as PrimaryHDU fast forward fixes to WCS fix forward merging of display --- package/FOC_reduction.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/package/FOC_reduction.py b/package/FOC_reduction.py index b42f57f..dcea9fa 100755 --- a/package/FOC_reduction.py +++ b/package/FOC_reduction.py @@ -273,7 +273,7 @@ def main(target=None, proposal_id=None, infiles=None, output_dir="./data", crop= flux_data=np.array([flux_data, flux_error, flux_mask]), flux_head=flux_head, ) - outfiles.append("/".join([data_folder, Stokes_hdul[0].header["FILENAME"] + ".fits"])) + outfiles.append("/".join([data_folder, Stokes_hdul["I_STOKES"].header["FILENAME"] + ".fits"])) # Step 5: # crop to desired region of interest (roi) From 427c4997a40bfdd784aa0745ed52236aa86f61c6 Mon Sep 17 00:00:00 2001 From: Thibault Barnouin Date: Wed, 19 Mar 2025 16:38:36 +0100 Subject: [PATCH 09/33] save raw flux in fits file and display forward fixes to display with WCS --- package/FOC_reduction.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/package/FOC_reduction.py b/package/FOC_reduction.py index dcea9fa..899c9e6 100755 --- a/package/FOC_reduction.py +++ b/package/FOC_reduction.py @@ -20,6 +20,7 @@ import lib.reduction as proj_red # Functions used in reduction pipeline import numpy as np from lib.utils import princ_angle, sci_not from matplotlib.colors import LogNorm +from astropy.wcs import WCS def main(target=None, proposal_id=None, infiles=None, output_dir="./data", crop=False, interactive=False): @@ -273,7 +274,7 @@ def main(target=None, proposal_id=None, infiles=None, output_dir="./data", crop= flux_data=np.array([flux_data, flux_error, flux_mask]), flux_head=flux_head, ) - outfiles.append("/".join([data_folder, Stokes_hdul["I_STOKES"].header["FILENAME"] + ".fits"])) + outfiles.append("/".join([data_folder, Stokes_hdul[0].header["FILENAME"] + ".fits"])) # Step 5: # crop to desired region of interest (roi) From e40f6bcf64a2e2c0148f8d8318f0a0ebc083ab11 Mon Sep 17 00:00:00 2001 From: Thibault Barnouin Date: Mon, 7 Apr 2025 17:20:00 +0200 Subject: [PATCH 10/33] modifications to add raw image display --- package/FOC_reduction.py | 10 +++++----- package/__init__.py | 2 +- package/lib/background.py | 3 ++- package/lib/plots.py | 24 +++++++++++++----------- 4 files changed, 21 insertions(+), 18 deletions(-) diff --git a/package/FOC_reduction.py b/package/FOC_reduction.py index 899c9e6..314e2ff 100755 --- a/package/FOC_reduction.py +++ b/package/FOC_reduction.py @@ -46,8 +46,8 @@ def main(target=None, proposal_id=None, infiles=None, output_dir="./data", crop= display_bkg = False # Data binning - pxsize = 4 - pxscale = "px" # pixel, arcsec or full + pxsize = 0.10 + pxscale = "arcsec" # pixel, arcsec or full rebin_operation = "sum" # sum or average # Alignement @@ -60,8 +60,8 @@ def main(target=None, proposal_id=None, infiles=None, output_dir="./data", crop= # Smoothing smoothing_function = "combine" # gaussian_after, weighted_gaussian_after, gaussian, weighted_gaussian or combine - smoothing_FWHM = 1.5 # If None, no smoothing is done - smoothing_scale = "px" # pixel or arcsec + smoothing_FWHM = 0.150 # If None, no smoothing is done + smoothing_scale = "arcsec" # pixel or arcsec # Rotation rotate_data = False @@ -71,7 +71,7 @@ def main(target=None, proposal_id=None, infiles=None, output_dir="./data", crop= P_cut = 3 # if >=1.0 cut on the signal-to-noise else cut on the confidence level in Q, U SNRi_cut = 1.0 # I measurments with SNR>30, which implies an uncertainty in P of 4.7%. flux_lim = None # lowest and highest flux displayed on plot, defaults to bkg and maximum in cut if None - scale_vec = 3 + scale_vec = 5 step_vec = 1 # plot all vectors in the array. if step_vec = 2, then every other vector will be plotted if step_vec = 0 then all vectors are displayed at full length # Pipeline start diff --git a/package/__init__.py b/package/__init__.py index a8ca1bd..59ef14c 100644 --- a/package/__init__.py +++ b/package/__init__.py @@ -1,2 +1,2 @@ from .lib import * -from . import FOC_reduction +from .FOC_reduction import main diff --git a/package/lib/background.py b/package/lib/background.py index 96509d2..3e8405f 100755 --- a/package/lib/background.py +++ b/package/lib/background.py @@ -18,7 +18,6 @@ import matplotlib.dates as mdates import matplotlib.pyplot as plt import numpy as np from astropy.time import Time -from lib.plots import plot_obs from matplotlib.colors import LogNorm from matplotlib.patches import Rectangle from scipy.optimize import curve_fit @@ -136,6 +135,8 @@ def display_bkg(data, background, std_bkg, headers, histograms=None, binning=Non fig2.subplots_adjust(hspace=0, wspace=0, right=1.0) fig2.colorbar(im2, ax=ax2, location="right", shrink=0.60, aspect=50, pad=0.025, label=r"Flux [$ergs \cdot cm^{-2} \cdot s^{-1} \cdot \AA^{-1}$]") + from .plots import plot_obs + if savename is not None: this_savename = deepcopy(savename) if savename[-4:] not in [".png", ".jpg", ".pdf"]: diff --git a/package/lib/plots.py b/package/lib/plots.py index a325b51..3797cbb 100755 --- a/package/lib/plots.py +++ b/package/lib/plots.py @@ -533,7 +533,7 @@ def polarization_map( cmap=kwargs["cmap"], alpha=1.0, ) - fig.colorbar(im, ax=ax, aspect=50, shrink=0.60, pad=0.025, label=r"$F_{\lambda}$ [$ergs \cdot cm^{-2} \cdot s^{-1} \cdot \AA$]") + fig.colorbar(im, ax=ax, aspect=50, shrink=0.60, pad=0.015, fraction=0.03, label=r"$F_{\lambda}$ [$ergs \cdot cm^{-2} \cdot s^{-1} \cdot \AA$]") levelsF = np.array([0.8, 2.0, 5.0, 10.0, 20.0, 50.0]) / 100.0 * vmax print("Flux density contour levels : ", levelsF) ax.contour(flux, levels=levelsF, transform=ax.get_transform(WCS(Stokes[0].header).celestial), colors="grey", linewidths=0.5) @@ -550,7 +550,7 @@ def polarization_map( vmin, vmax = 1.0 / 2.0 * np.median(np.sqrt(stk_cov[0, 0][stkI > 0.0]) * convert_flux), np.max(stkI[stkI > 0.0] * convert_flux) im = ax.imshow(stkI * convert_flux, norm=LogNorm(vmin, vmax), aspect="equal", cmap=kwargs["cmap"], alpha=1.0) ax.plot(*WCS(Stokes[1]).wcs.crpix, "g+") - fig.colorbar(im, ax=ax, aspect=50, shrink=0.60, pad=0.025, label=r"$F_{\lambda}$ [$ergs \cdot cm^{-2} \cdot s^{-1} \cdot \AA^{-1}$]") + fig.colorbar(im, ax=ax, aspect=50, shrink=0.60, pad=0.015, fraction=0.03, label=r"$F_{\lambda}$ [$ergs \cdot cm^{-2} \cdot s^{-1} \cdot \AA^{-1}$]") # levelsI = np.array([0.8, 2.0, 5.0, 10.0, 20.0, 50.0]) / 100.0 * vmax # print("Stokes I contour levels : ", levelsI) # ax.contour(stkI * convert_flux, levels=levelsI, colors="grey", linewidths=0.5) @@ -569,7 +569,9 @@ def polarization_map( vmin, vmax = 1.0 / 2.0 * np.median(np.sqrt(stk_cov[0, 0][stkI > 0.0]) * convert_flux), np.max(stkI[stkI > 0.0] * convert_flux) pfmax = (stkI[stkI > 0.0] * pol[stkI > 0.0] * convert_flux).max() im = ax.imshow(stkI * convert_flux * pol, norm=LogNorm(vmin, vmax), aspect="equal", cmap=kwargs["cmap"], alpha=1.0) - fig.colorbar(im, ax=ax, aspect=50, shrink=0.60, pad=0.025, label=r"$F_{\lambda} \cdot P$ [$ergs \cdot cm^{-2} \cdot s^{-1} \cdot \AA^{-1}$]") + fig.colorbar( + im, ax=ax, aspect=50, shrink=0.60, pad=0.015, fraction=0.03, label=r"$F_{\lambda} \cdot P$ [$ergs \cdot cm^{-2} \cdot s^{-1} \cdot \AA^{-1}$]" + ) # levelsPf = np.linspace(0.0.60, 0.50, 5) * pfmax levelsPf = np.array([1.73, 13.0, 33.0, 66.0]) / 100.0 * pfmax print("Polarized flux density contour levels : ", levelsPf) @@ -579,13 +581,13 @@ def polarization_map( display = "p" vmin, vmax = 0.0, 100.0 im = ax.imshow(pol * 100.0, vmin=vmin, vmax=vmax, aspect="equal", cmap=kwargs["cmap"], alpha=1.0) - fig.colorbar(im, ax=ax, aspect=50, shrink=0.60, pad=0.025, label=r"$P$ [%]") + fig.colorbar(im, ax=ax, aspect=50, shrink=0.60, pad=0.015, fraction=0.03, label=r"$P$ [%]") elif display.lower() in ["pa", "pang", "pol_ang"]: # Display polarization degree map display = "pa" vmin, vmax = 0.0, 180.0 im = ax.imshow(princ_angle(pang), vmin=vmin, vmax=vmax, aspect="equal", cmap=kwargs["cmap"], alpha=1.0) - fig.colorbar(im, ax=ax, aspect=50, shrink=0.60, pad=0.025, label=r"$\theta_P$ [°]") + fig.colorbar(im, ax=ax, aspect=50, shrink=0.60, pad=0.015, fraction=0.03, label=r"$\theta_P$ [°]") elif display.lower() in ["s_p", "pol_err", "pol_deg_err"]: # Display polarization degree error map display = "s_p" @@ -595,7 +597,7 @@ def polarization_map( else: vmin, vmax = 0.0, 100.0 im = ax.imshow(pol_err * 100.0, vmin=vmin, vmax=vmax, aspect="equal", cmap="inferno_r", alpha=1.0) - fig.colorbar(im, ax=ax, aspect=50, shrink=0.60, pad=0.025, label=r"$\sigma_P$ [%]") + fig.colorbar(im, ax=ax, aspect=50, shrink=0.60, pad=0.015, fraction=0.03, label=r"$\sigma_P$ [%]") elif display.lower() in ["s_i", "i_err"]: # Display intensity error map display = "s_i" @@ -607,7 +609,7 @@ def polarization_map( im = ax.imshow(np.sqrt(stk_cov[0, 0]) * convert_flux, norm=LogNorm(vmin, vmax), aspect="equal", cmap="inferno_r", alpha=1.0) else: im = ax.imshow(np.sqrt(stk_cov[0, 0]) * convert_flux, aspect="equal", cmap=kwargs["cmap"], alpha=1.0) - fig.colorbar(im, ax=ax, aspect=50, shrink=0.60, pad=0.025, label=r"$\sigma_I$ [$ergs \cdot cm^{-2} \cdot s^{-1} \cdot \AA^{-1}$]") + fig.colorbar(im, ax=ax, aspect=50, shrink=0.60, pad=0.015, fraction=0.03, label=r"$\sigma_I$ [$ergs \cdot cm^{-2} \cdot s^{-1} \cdot \AA^{-1}$]") elif display.lower() in ["snri"]: # Display I_stokes signal-to-noise map display = "snri" @@ -619,7 +621,7 @@ def polarization_map( ax.contour(SNRi, levels=levelsSNRi, colors="grey", linewidths=0.5) else: im = ax.imshow(SNRi, aspect="equal", cmap=kwargs["cmap"], alpha=1.0) - fig.colorbar(im, ax=ax, aspect=50, shrink=0.60, pad=0.025, label=r"$I_{Stokes}/\sigma_{I}$") + fig.colorbar(im, ax=ax, aspect=50, shrink=0.60, pad=0.015, fraction=0.03, label=r"$I_{Stokes}/\sigma_{I}$") elif display.lower() in ["snr", "snrp"]: # Display polarization degree signal-to-noise map display = "snrp" @@ -631,7 +633,7 @@ def polarization_map( ax.contour(SNRp, levels=levelsSNRp, colors="grey", linewidths=0.5) else: im = ax.imshow(SNRp, aspect="equal", cmap=kwargs["cmap"], alpha=1.0) - fig.colorbar(im, ax=ax, aspect=50, shrink=0.60, pad=0.025, label=r"$P/\sigma_{P}$") + fig.colorbar(im, ax=ax, aspect=50, shrink=0.60, pad=0.015, fraction=0.03, label=r"$P/\sigma_{P}$") elif display.lower() in ["conf", "confp"]: # Display polarization degree signal-to-noise map display = "confp" @@ -640,7 +642,7 @@ def polarization_map( levelsconfp = np.array([0.500, 0.900, 0.990, 0.999]) print("confp contour levels : ", levelsconfp) ax.contour(confP, levels=levelsconfp, colors="grey", linewidths=0.5) - fig.colorbar(im, ax=ax, aspect=50, shrink=0.60, pad=0.025, label=r"$Conf_{P}$") + fig.colorbar(im, ax=ax, aspect=50, shrink=0.60, pad=0.015, fraction=0.03, label=r"$Conf_{P}$") else: # Defaults to intensity map if flux_lim is not None: @@ -655,7 +657,7 @@ def polarization_map( cmap=kwargs["cmap"], alpha=1.0, ) - fig.colorbar(im, ax=ax, aspect=50, shrink=0.60, pad=0.025, label=r"$F_{\lambda}$ [$ergs \cdot cm^{-2} \cdot s^{-1} \cdot \AA$]") + fig.colorbar(im, ax=ax, aspect=50, shrink=0.60, pad=0.015, fraction=0.03, label=r"$F_{\lambda}$ [$ergs \cdot cm^{-2} \cdot s^{-1} \cdot \AA$]") I_diluted, I_diluted_err = np.sum(flux[flux_mask]), np.sqrt(np.sum(flux_error[flux_mask] ** 2)) plt.rcParams.update({"font.size": 11}) From 40982f1fe07d2a1ba9da0d1f8d70fb4e9db41582 Mon Sep 17 00:00:00 2001 From: Thibault Barnouin Date: Tue, 11 Mar 2025 16:14:14 +0100 Subject: [PATCH 11/33] small improvments to ConfCenter et emission_center --- package/lib/utils.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/package/lib/utils.py b/package/lib/utils.py index 9e36dbc..e56a29e 100755 --- a/package/lib/utils.py +++ b/package/lib/utils.py @@ -109,7 +109,7 @@ def CenterConf(mask, PA, sPA): conf[np.isfinite(PA)] = 0.5 * np.exp(-0.5 * chi2[np.isfinite(PA)]) c0 = np.unravel_index(np.argmax(conf), conf.shape)[::-1] - result = minimize(chisq, c0, bounds=[(0, PA.shape[1]), (0.0, PA.shape[0])], method="trust-constr") + result = minimize(chisq, c0, bounds=[(0, PA.shape[1]), (0.0, PA.shape[0])], method="TNC") if result.success: print("Center of emission found: reduced chi_squared {0:.2f}/{1:d}={2:.2f}".format(chisq(result.x) * (Nobs - 2), Nobs - 2, chisq(result.x))) else: From 78eb00c4867afba7f58f1537bd2d9780ebbc77dc Mon Sep 17 00:00:00 2001 From: Thibault Barnouin Date: Fri, 14 Mar 2025 14:20:38 +0100 Subject: [PATCH 12/33] rollback CenterConf to simple CDF of chi2 for 2 dof --- package/lib/utils.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/package/lib/utils.py b/package/lib/utils.py index e56a29e..9e36dbc 100755 --- a/package/lib/utils.py +++ b/package/lib/utils.py @@ -109,7 +109,7 @@ def CenterConf(mask, PA, sPA): conf[np.isfinite(PA)] = 0.5 * np.exp(-0.5 * chi2[np.isfinite(PA)]) c0 = np.unravel_index(np.argmax(conf), conf.shape)[::-1] - result = minimize(chisq, c0, bounds=[(0, PA.shape[1]), (0.0, PA.shape[0])], method="TNC") + result = minimize(chisq, c0, bounds=[(0, PA.shape[1]), (0.0, PA.shape[0])], method="trust-constr") if result.success: print("Center of emission found: reduced chi_squared {0:.2f}/{1:d}={2:.2f}".format(chisq(result.x) * (Nobs - 2), Nobs - 2, chisq(result.x))) else: From 9be8e28cc97147bc6e1869a2f977593f1b03f9cb Mon Sep 17 00:00:00 2001 From: Thibault Barnouin Date: Tue, 1 Apr 2025 17:02:31 +0200 Subject: [PATCH 13/33] fix WCS computation, cdelt should not be sorted fix rebase display on main --- package/FOC_reduction.py | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/package/FOC_reduction.py b/package/FOC_reduction.py index 314e2ff..d3aee2c 100755 --- a/package/FOC_reduction.py +++ b/package/FOC_reduction.py @@ -46,8 +46,8 @@ def main(target=None, proposal_id=None, infiles=None, output_dir="./data", crop= display_bkg = False # Data binning - pxsize = 0.10 - pxscale = "arcsec" # pixel, arcsec or full + pxsize = 4 + pxscale = "px" # pixel, arcsec or full rebin_operation = "sum" # sum or average # Alignement @@ -60,8 +60,8 @@ def main(target=None, proposal_id=None, infiles=None, output_dir="./data", crop= # Smoothing smoothing_function = "combine" # gaussian_after, weighted_gaussian_after, gaussian, weighted_gaussian or combine - smoothing_FWHM = 0.150 # If None, no smoothing is done - smoothing_scale = "arcsec" # pixel or arcsec + smoothing_FWHM = 1.5 # If None, no smoothing is done + smoothing_scale = "px" # pixel or arcsec # Rotation rotate_data = False From cdc19c42b4d999ff8ada82ffc660143b008a9972 Mon Sep 17 00:00:00 2001 From: Thibault Barnouin Date: Mon, 7 Apr 2025 17:02:32 +0200 Subject: [PATCH 14/33] fix error on WCS rotation that could mirror the image --- package/lib/fits.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/package/lib/fits.py b/package/lib/fits.py index 4ea9a47..67a62d9 100755 --- a/package/lib/fits.py +++ b/package/lib/fits.py @@ -16,7 +16,7 @@ from astropy.io import fits from astropy.wcs import WCS from .convex_hull import clean_ROI -from .utils import wcs_PA, princ_angle +from .utils import wcs_PA def get_obs_data(infiles, data_folder="", compute_flux=False): @@ -78,7 +78,7 @@ def get_obs_data(infiles, data_folder="", compute_flux=False): for key, val in new_wcs.to_header().items(): header[key] = val try: - header["ORIENTAT"] = princ_angle(float(header["ORIENTAT"])) + header["ORIENTAT"] = float(header["ORIENTAT"]) except KeyError: header["ORIENTAT"] = wcs_PA(new_wcs.wcs.pc[1, 0], np.diag(new_wcs.wcs.pc).mean()) From e9efd02e3fce35194f4975c5eaef2bceebc8fec9 Mon Sep 17 00:00:00 2001 From: Thibault Barnouin Date: Fri, 14 Mar 2025 14:30:30 +0100 Subject: [PATCH 15/33] Save the raw total flux image as PrimaryHDU fix for rebase display on main fix rebase display on main --- package/FOC_reduction.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/package/FOC_reduction.py b/package/FOC_reduction.py index d3aee2c..5abef77 100755 --- a/package/FOC_reduction.py +++ b/package/FOC_reduction.py @@ -274,7 +274,7 @@ def main(target=None, proposal_id=None, infiles=None, output_dir="./data", crop= flux_data=np.array([flux_data, flux_error, flux_mask]), flux_head=flux_head, ) - outfiles.append("/".join([data_folder, Stokes_hdul[0].header["FILENAME"] + ".fits"])) + outfiles.append("/".join([data_folder, Stokes_hdul["I_STOKES"].header["FILENAME"] + ".fits"])) # Step 5: # crop to desired region of interest (roi) From 05f800b2820c55d3fe352e87eea6837a1743303b Mon Sep 17 00:00:00 2001 From: Thibault Barnouin Date: Wed, 19 Mar 2025 16:38:36 +0100 Subject: [PATCH 16/33] save raw flux in fits file and display fix rebase display on main fix rebase display on main --- package/FOC_reduction.py | 2 +- package/lib/plots.py | 4 ++-- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/package/FOC_reduction.py b/package/FOC_reduction.py index 5abef77..d3aee2c 100755 --- a/package/FOC_reduction.py +++ b/package/FOC_reduction.py @@ -274,7 +274,7 @@ def main(target=None, proposal_id=None, infiles=None, output_dir="./data", crop= flux_data=np.array([flux_data, flux_error, flux_mask]), flux_head=flux_head, ) - outfiles.append("/".join([data_folder, Stokes_hdul["I_STOKES"].header["FILENAME"] + ".fits"])) + outfiles.append("/".join([data_folder, Stokes_hdul[0].header["FILENAME"] + ".fits"])) # Step 5: # crop to desired region of interest (roi) diff --git a/package/lib/plots.py b/package/lib/plots.py index 3797cbb..994b6b6 100755 --- a/package/lib/plots.py +++ b/package/lib/plots.py @@ -3575,7 +3575,7 @@ class pol_map(object): @property def wcs(self): - return WCS(self.Stokes["I_STOKES"].header).celestial.deepcopy() + return WCS(self.Stokes[0].header).celestial.deepcopy() @property def Flux(self): @@ -3735,7 +3735,7 @@ class pol_map(object): back_length=0.0, head_length=10.0, head_width=10.0, - angle=-self.Stokes["I_STOKES"].header["orientat"], + angle=-self.Stokes[0].header["orientat"], color="white", text_props={"ec": None, "fc": "w", "alpha": 1, "lw": 0.4}, arrow_props={"ec": None, "fc": "w", "alpha": 1, "lw": 1}, From 29786c2fa4957bab3fd097e3a10b390818bc317f Mon Sep 17 00:00:00 2001 From: Thibault Barnouin Date: Tue, 8 Apr 2025 19:48:46 +0200 Subject: [PATCH 17/33] some more formatting --- package/FOC_reduction.py | 7 +- package/lib/plots.py | 996 ++++++--------------------------------- package/lib/reduction.py | 31 +- 3 files changed, 141 insertions(+), 893 deletions(-) diff --git a/package/FOC_reduction.py b/package/FOC_reduction.py index d3aee2c..d3782cb 100755 --- a/package/FOC_reduction.py +++ b/package/FOC_reduction.py @@ -290,12 +290,7 @@ def main(target=None, proposal_id=None, infiles=None, output_dir="./data", crop= print( "F_int({0:.0f} Angs) = ({1} ± {2})e{3} ergs.cm^-2.s^-1.Angs^-1".format( flux_head["PHOTPLAM"], - *sci_not( - flux_data[flux_mask].sum() * flux_head["PHOTFLAM"], - np.sqrt(np.sum(flux_error[flux_mask] ** 2)) * flux_head["PHOTFLAM"], - 2, - out=int, - ), + *sci_not(flux_data[flux_mask].sum() * flux_head["PHOTFLAM"], np.sqrt(np.sum(flux_error[flux_mask] ** 2)) * flux_head["PHOTFLAM"], 2, out=int), ) ) print("P_int = {0:.1f} ± {1:.1f} %".format(header_stokes["p_int"] * 100.0, np.ceil(header_stokes["sP_int"] * 1000.0) / 10.0)) diff --git a/package/lib/plots.py b/package/lib/plots.py index 994b6b6..931abe6 100755 --- a/package/lib/plots.py +++ b/package/lib/plots.py @@ -54,10 +54,7 @@ from matplotlib.colors import LogNorm from matplotlib.patches import Circle, FancyArrowPatch, Rectangle from matplotlib.path import Path from matplotlib.widgets import Button, LassoSelector, RectangleSelector, Slider, TextBox -from mpl_toolkits.axes_grid1.anchored_artists import ( - AnchoredDirectionArrows, - AnchoredSizeBar, -) +from mpl_toolkits.axes_grid1.anchored_artists import AnchoredDirectionArrows, AnchoredSizeBar try: from .utils import PCconf, princ_angle, rot2D, sci_not @@ -65,15 +62,7 @@ except ImportError: from utils import PCconf, princ_angle, rot2D, sci_not -def plot_obs( - data_array, - headers, - rectangle=None, - shifts=None, - savename=None, - plots_folder="", - **kwargs, -): +def plot_obs(data_array, headers, rectangle=None, shifts=None, savename=None, plots_folder="", **kwargs): """ Plots raw observation imagery with some information on the instrument and filters. @@ -104,14 +93,7 @@ def plot_obs( plt.rcParams.update({"font.size": 10}) nb_obs = np.max([np.sum([head["filtnam1"] == curr_pol for head in headers]) for curr_pol in ["POL0", "POL60", "POL120"]]) shape = np.array((3, nb_obs)) - fig, ax = plt.subplots( - shape[0], - shape[1], - figsize=(3 * shape[1], 3 * shape[0]), - layout="compressed", - sharex=True, - sharey=True, - ) + fig, ax = plt.subplots(shape[0], shape[1], figsize=(3 * shape[1], 3 * shape[0]), layout="compressed", sharex=True, sharey=True) used = np.zeros(shape, dtype=bool) r_pol = dict(pol0=0, pol60=1, pol120=2) c_pol = dict(pol0=0, pol60=0, pol120=0) @@ -134,14 +116,8 @@ def plot_obs( vmin, vmax = kwargs["vmin"], kwargs["vmax"] del kwargs["vmin"], kwargs["vmax"] else: - vmin, vmax = ( - convert * data[data > 0.0].min() / 10.0, - convert * data[data > 0.0].max(), - ) - for key, value in [ - ["cmap", [["cmap", "inferno"]]], - ["norm", [["norm", LogNorm(vmin, vmax)]]], - ]: + vmin, vmax = (convert * data[data > 0.0].min() / 10.0, convert * data[data > 0.0].max()) + for key, value in [["cmap", [["cmap", "inferno"]]], ["norm", [["norm", LogNorm(vmin, vmax)]]]]: try: _ = kwargs[key] except KeyError: @@ -153,35 +129,12 @@ def plot_obs( if shifts is not None: x, y = np.array(data.shape[::-1]) / 2.0 - shifts[i] dx, dy = shifts[i] - ax_curr.arrow( - x, - y, - dx, - dy, - length_includes_head=True, - width=0.1, - head_width=0.3, - color="g", - ) + ax_curr.arrow(x, y, dx, dy, length_includes_head=True, width=0.1, head_width=0.3, color="g") ax_curr.plot([x, x], [0, data.shape[0] - 1], "--", lw=2, color="g", alpha=0.85) ax_curr.plot([0, data.shape[1] - 1], [y, y], "--", lw=2, color="g", alpha=0.85) # position of centroid - ax_curr.plot( - [data.shape[1] / 2, data.shape[1] / 2], - [0, data.shape[0] - 1], - "--", - lw=2, - color="b", - alpha=0.85, - ) - ax_curr.plot( - [0, data.shape[1] - 1], - [data.shape[0] / 2, data.shape[0] / 2], - "--", - lw=2, - color="b", - alpha=0.85, - ) + ax_curr.plot([data.shape[1] / 2, data.shape[1] / 2], [0, data.shape[0] - 1], "--", lw=2, color="b", alpha=0.85) + ax_curr.plot([0, data.shape[1] - 1], [data.shape[0] / 2, data.shape[0] / 2], "--", lw=2, color="b", alpha=0.85) cr_x, cr_y = head["CRPIX1"], head["CRPIX2"] # Plot WCS reference point ax_curr.plot([cr_x], [cr_y], "+", lw=2, color="r", alpha=0.85) @@ -189,31 +142,11 @@ def plot_obs( x, y, width, height, angle, color = rectangle[i] ax_curr.add_patch(Rectangle((x, y), width, height, angle=angle, edgecolor=color, fill=False)) ax_curr.annotate( - instr + ":" + rootname, - color="white", - fontsize=10, - xy=(0.01, 1.00), - xycoords="axes fraction", - verticalalignment="top", - horizontalalignment="left", + instr + ":" + rootname, color="white", fontsize=10, xy=(0.01, 1.00), xycoords="axes fraction", verticalalignment="top", horizontalalignment="left" ) + ax_curr.annotate(filt, color="white", fontsize=15, xy=(0.01, 0.01), xycoords="axes fraction", verticalalignment="bottom", horizontalalignment="left") ax_curr.annotate( - filt, - color="white", - fontsize=15, - xy=(0.01, 0.01), - xycoords="axes fraction", - verticalalignment="bottom", - horizontalalignment="left", - ) - ax_curr.annotate( - exptime, - color="white", - fontsize=10, - xy=(1.00, 0.01), - xycoords="axes fraction", - verticalalignment="bottom", - horizontalalignment="right", + exptime, color="white", fontsize=10, xy=(1.00, 0.01), xycoords="axes fraction", verticalalignment="bottom", horizontalalignment="right" ) unused = np.logical_not(used) ii, jj = np.indices(shape) @@ -221,26 +154,13 @@ def plot_obs( fig.delaxes(ax[i][j]) # fig.subplots_adjust(hspace=0.01, wspace=0.01, right=1.02) - fig.colorbar( - im, - ax=ax, - location="right", - shrink=0.75, - aspect=50, - pad=0.025, - label=r"$F_{\lambda}$ [$ergs \cdot cm^{-2} \cdot s^{-1} \cdot \AA^{-1}$]", - ) + fig.colorbar(im, ax=ax, location="right", shrink=0.75, aspect=50, pad=0.025, label=r"$F_{\lambda}$ [$ergs \cdot cm^{-2} \cdot s^{-1} \cdot \AA^{-1}$]") if savename is not None: # fig.suptitle(savename) if savename[-4:] not in [".png", ".jpg", ".pdf"]: savename += ".pdf" - fig.savefig( - path_join(plots_folder, savename), - bbox_inches="tight", - dpi=150, - facecolor="None", - ) + fig.savefig(path_join(plots_folder, savename), bbox_inches="tight", dpi=150, facecolor="None") plt.show() return 0 @@ -299,12 +219,7 @@ def plot_Stokes(Stokes, savename=None, plots_folder=""): savename += "_IQU.pdf" else: savename = savename[:-4] + "_IQU" + savename[-4:] - fig.savefig( - path_join(plots_folder, savename), - bbox_inches="tight", - dpi=150, - facecolor="None", - ) + fig.savefig(path_join(plots_folder, savename), bbox_inches="tight", dpi=150, facecolor="None") return 0 @@ -408,10 +323,7 @@ def polarization_map( # Compute confidence level map QN, UN, QN_ERR, UN_ERR = np.full((4, stkI.shape[0], stkI.shape[1]), np.nan) - for nflux, sflux in zip( - [QN, UN, QN_ERR, UN_ERR], - [stkQ, stkU, np.sqrt(stk_cov[1, 1]), np.sqrt(stk_cov[2, 2])], - ): + for nflux, sflux in zip([QN, UN, QN_ERR, UN_ERR], [stkQ, stkU, np.sqrt(stk_cov[1, 1]), np.sqrt(stk_cov[2, 2])]): nflux[stkI > 0.0] = sflux[stkI > 0.0] / stkI[stkI > 0.0] confP = PCconf(QN, UN, QN_ERR, UN_ERR) @@ -476,11 +388,7 @@ def polarization_map( ax.set_ylabel("Declination (J2000)", labelpad=-1) vmin, vmax = 0.0, stkI.max() * convert_flux - for key, value in [ - ["cmap", [["cmap", "inferno"]]], - ["width", [["width", 0.4]]], - ["linewidth", [["linewidth", 0.6]]], - ]: + for key, value in [["cmap", [["cmap", "inferno"]]], ["width", [["width", 0.4]]], ["linewidth", [["linewidth", 0.6]]]]: try: _ = kwargs[key] except KeyError: @@ -526,12 +434,7 @@ def polarization_map( imflux, cr = flux.copy(), WCS(Stokes[0].header).wcs.crpix.astype(int) imflux[cr[1] - 1 : cr[1] + 2, cr[0] - 1 : cr[0] + 2] = np.nan im = ax.imshow( - imflux, - transform=ax.get_transform(WCS(Stokes[0].header).celestial), - norm=LogNorm(vmin, vmax), - aspect="equal", - cmap=kwargs["cmap"], - alpha=1.0, + imflux, transform=ax.get_transform(WCS(Stokes[0].header).celestial), norm=LogNorm(vmin, vmax), aspect="equal", cmap=kwargs["cmap"], alpha=1.0 ) fig.colorbar(im, ax=ax, aspect=50, shrink=0.60, pad=0.015, fraction=0.03, label=r"$F_{\lambda}$ [$ergs \cdot cm^{-2} \cdot s^{-1} \cdot \AA$]") levelsF = np.array([0.8, 2.0, 5.0, 10.0, 20.0, 50.0]) / 100.0 * vmax @@ -662,18 +565,7 @@ def polarization_map( plt.rcParams.update({"font.size": 11}) px_size = wcs.wcs.get_cdelt()[0] * 3600.0 - px_sc = AnchoredSizeBar( - ax.transData, - 1.0 / px_size, - "1 arcsec", - 3, - pad=0.25, - sep=5, - borderpad=0.25, - frameon=False, - size_vertical=0.005, - color=font_color, - ) + px_sc = AnchoredSizeBar(ax.transData, 1.0 / px_size, "1 arcsec", 3, pad=0.25, sep=5, borderpad=0.25, frameon=False, size_vertical=0.005, color=font_color) north_dir = AnchoredDirectionArrows( ax.transAxes, "E", @@ -698,10 +590,7 @@ def polarization_map( step_vec = 1 scale_vec = 2.0 X, Y = np.meshgrid(np.arange(stkI.shape[1]), np.arange(stkI.shape[0])) - U, V = ( - poldata * np.cos(np.pi / 2.0 + pangdata * np.pi / 180.0), - poldata * np.sin(np.pi / 2.0 + pangdata * np.pi / 180.0), - ) + U, V = (poldata * np.cos(np.pi / 2.0 + pangdata * np.pi / 180.0), poldata * np.sin(np.pi / 2.0 + pangdata * np.pi / 180.0)) ax.quiver( X[::step_vec, ::step_vec], Y[::step_vec, ::step_vec], @@ -721,16 +610,7 @@ def polarization_map( edgecolor="k", ) pol_sc = AnchoredSizeBar( - ax.transData, - scale_vec, - r"$P$= 100 %", - 4, - pad=0.25, - sep=5, - borderpad=0.25, - frameon=False, - size_vertical=0.005, - color=font_color, + ax.transData, scale_vec, r"$P$= 100 %", 4, pad=0.25, sep=5, borderpad=0.25, frameon=False, size_vertical=0.005, color=font_color ) ax.add_artist(pol_sc) @@ -739,8 +619,7 @@ def polarization_map( ax.annotate( r"$F_{{\lambda}}^{{int}}$({0:.0f} $\AA$) = {1} $ergs \cdot cm^{{-2}} \cdot s^{{-1}} \cdot \AA^{{-1}}$".format( - pivot_wav, - sci_not(I_diluted, I_diluted_err, 2), + pivot_wav, sci_not(I_diluted, I_diluted_err, 2) ) + "\n" + r"$P^{{int}}$ = {0:.1f} $\pm$ {1:.1f} %".format(P_diluted * 100.0, P_diluted_err * 100.0) @@ -759,8 +638,7 @@ def polarization_map( ax.add_artist(north_dir) ax.annotate( r"$F_{{\lambda}}^{{int}}$({0:.0f} $\AA$) = {1} $ergs \cdot cm^{{-2}} \cdot s^{{-1}} \cdot \AA^{{-1}}$".format( - pivot_wav, - sci_not(I_diluted, I_diluted_err, 2), + pivot_wav, sci_not(I_diluted, I_diluted_err, 2) ), color="white", xy=(0.01, 1.00), @@ -779,12 +657,7 @@ def polarization_map( if savename is not None: if savename[-4:] not in [".png", ".jpg", ".pdf"]: savename += ".pdf" - fig.savefig( - path_join(plots_folder, savename), - bbox_inches="tight", - dpi=300, - facecolor="None", - ) + fig.savefig(path_join(plots_folder, savename), bbox_inches="tight", dpi=300, facecolor="None") return fig, ax @@ -820,26 +693,14 @@ class align_maps(object): self.other_data = self.other_data[0] self.map_convert, self.map_unit = ( - ( - float(self.map_header["photflam"]), - r"$ergs \cdot cm^{-2} \cdot s^{-1} \cdot \AA^{-1}$", - ) + (float(self.map_header["photflam"]), r"$ergs \cdot cm^{-2} \cdot s^{-1} \cdot \AA^{-1}$") if "PHOTFLAM" in list(self.map_header.keys()) - else ( - 1.0, - self.map_header["bunit"] if "BUNIT" in list(self.map_header.keys()) else "Arbitray Units", - ) + else (1.0, self.map_header["bunit"] if "BUNIT" in list(self.map_header.keys()) else "Arbitray Units") ) self.other_convert, self.other_unit = ( - ( - float(self.other_header["photflam"]), - r"$ergs \cdot cm^{-2} \cdot s^{-1} \cdot \AA^{-1}$", - ) + (float(self.other_header["photflam"]), r"$ergs \cdot cm^{-2} \cdot s^{-1} \cdot \AA^{-1}$") if "PHOTFLAM" in list(self.other_header.keys()) - else ( - 1.0, - self.other_header["bunit"] if "BUNIT" in list(self.other_header.keys()) else "Arbitray Units", - ) + else (1.0, self.other_header["bunit"] if "BUNIT" in list(self.other_header.keys()) else "Arbitray Units") ) self.map_observer = ( "/".join([self.map_header["telescop"], self.map_header["instrume"]]) if "INSTRUME" in list(self.map_header.keys()) else self.map_header["telescop"] @@ -858,14 +719,8 @@ class align_maps(object): # Plot the UV map other_kwargs = deepcopy(kwargs) - vmin, vmax = ( - self.map_data[self.map_data > 0.0].max() / 1e3 * self.map_convert, - self.map_data[self.map_data > 0.0].max() * self.map_convert, - ) - for key, value in [ - ["cmap", [["cmap", "inferno"]]], - ["norm", [["norm", LogNorm(vmin, vmax)]]], - ]: + vmin, vmax = (self.map_data[self.map_data > 0.0].max() / 1e3 * self.map_convert, self.map_data[self.map_data > 0.0].max() * self.map_convert) + for key, value in [["cmap", [["cmap", "inferno"]]], ["norm", [["norm", LogNorm(vmin, vmax)]]]]: try: _ = kwargs[key] except KeyError: @@ -957,10 +812,7 @@ class align_maps(object): self.other_data[self.other_data > 0.0].max() / 1e3 * self.other_convert, self.other_data[self.other_data > 0.0].max() * self.other_convert, ) - for key, value in [ - ["cmap", [["cmap", "inferno"]]], - ["norm", [["norm", LogNorm(vmin, vmax)]]], - ]: + for key, value in [["cmap", [["cmap", "inferno"]]], ["norm", [["norm", LogNorm(vmin, vmax)]]]]: try: _ = other_kwargs[key] except KeyError: @@ -1064,10 +916,7 @@ class align_maps(object): else: self.map_wcs.wcs.crpix = np.array(self.cr_map.get_data()) + (1.0, 1.0) if np.array(self.cr_other.get_data()).shape == (2, 1): - self.other_wcs.wcs.crpix = np.array(self.cr_other.get_data())[:, 0] + ( - 1.0, - 1.0, - ) + self.other_wcs.wcs.crpix = np.array(self.cr_other.get_data())[:, 0] + (1.0, 1.0) else: self.other_wcs.wcs.crpix = np.array(self.cr_other.get_data()) + (1.0, 1.0) self.map_wcs.wcs.crval = np.array(self.map_wcs.pixel_to_world_values(*self.map_wcs.wcs.crpix)) @@ -1119,15 +968,7 @@ class overplot_radio(align_maps): Inherit from class align_maps in order to get the same WCS on both maps. """ - def overplot( - self, - levels=None, - P_cut=0.99, - SNRi_cut=1.0, - scale_vec=2, - savename=None, - **kwargs, - ): + def overplot(self, levels=None, P_cut=0.99, SNRi_cut=1.0, scale_vec=2, savename=None, **kwargs): self.Stokes_UV = self.map self.wcs_UV = self.map_wcs # Get Data @@ -1141,10 +982,7 @@ class overplot_radio(align_maps): pang = self.Stokes_UV["POL_ANG"].data # Compute confidence level map QN, UN, QN_ERR, UN_ERR = np.full((4, stkI.shape[0], stkI.shape[1]), np.nan) - for nflux, sflux in zip( - [QN, UN, QN_ERR, UN_ERR], - [stkQ, stkU, np.sqrt(stk_cov[1, 1]), np.sqrt(stk_cov[2, 2])], - ): + for nflux, sflux in zip([QN, UN, QN_ERR, UN_ERR], [stkQ, stkU, np.sqrt(stk_cov[1, 1]), np.sqrt(stk_cov[2, 2])]): nflux[stkI > 0.0] = sflux[stkI > 0.0] / stkI[stkI > 0.0] confP = PCconf(QN, UN, QN_ERR, UN_ERR) @@ -1176,14 +1014,8 @@ class overplot_radio(align_maps): self.fig_overplot.subplots_adjust(hspace=0, wspace=0, bottom=0.1, left=0.1, top=0.8, right=1) # Display UV intensity map with polarization vectors - vmin, vmax = ( - stkI[np.isfinite(stkI)].max() / 1e3 * self.map_convert, - stkI[np.isfinite(stkI)].max() * self.map_convert, - ) - for key, value in [ - ["cmap", [["cmap", "inferno"]]], - ["norm", [["norm", LogNorm(vmin, vmax)]]], - ]: + vmin, vmax = (stkI[np.isfinite(stkI)].max() / 1e3 * self.map_convert, stkI[np.isfinite(stkI)].max() * self.map_convert) + for key, value in [["cmap", [["cmap", "inferno"]]], ["norm", [["norm", LogNorm(vmin, vmax)]]]]: try: _ = kwargs[key] except KeyError: @@ -1218,19 +1050,9 @@ class overplot_radio(align_maps): else: self.ax_overplot.set_facecolor("white") font_color = "black" - self.im = self.ax_overplot.imshow( - stkI * self.map_convert, - aspect="equal", - label="{0:s} observation".format(self.map_observer), - **kwargs, - ) + self.im = self.ax_overplot.imshow(stkI * self.map_convert, aspect="equal", label="{0:s} observation".format(self.map_observer), **kwargs) self.cbar = self.fig_overplot.colorbar( - self.im, - ax=self.ax_overplot, - aspect=50, - shrink=0.75, - pad=0.025, - label=r"$F_{{\lambda}}$ [{0:s}]".format(self.map_unit), + self.im, ax=self.ax_overplot, aspect=50, shrink=0.75, pad=0.025, label=r"$F_{{\lambda}}$ [{0:s}]".format(self.map_unit) ) # Display full size polarization vectors @@ -1241,10 +1063,7 @@ class overplot_radio(align_maps): self.scale_vec = scale_vec step_vec = 1 self.X, self.Y = np.meshgrid(np.arange(stkI.shape[1]), np.arange(stkI.shape[0])) - self.U, self.V = ( - pol * np.cos(np.pi / 2.0 + pang * np.pi / 180.0), - pol * np.sin(np.pi / 2.0 + pang * np.pi / 180.0), - ) + self.U, self.V = (pol * np.cos(np.pi / 2.0 + pang * np.pi / 180.0), pol * np.sin(np.pi / 2.0 + pang * np.pi / 180.0)) self.Q = self.ax_overplot.quiver( self.X[::step_vec, ::step_vec], self.Y[::step_vec, ::step_vec], @@ -1281,11 +1100,7 @@ class overplot_radio(align_maps): self.ax_overplot.set_ylabel(label="Declination (J2000)", labelpad=-1) self.fig_overplot.suptitle( "{0:s} polarization map of {1:s} overplotted with {2:s} {3:.2f}GHz map in {4:s}.".format( - self.map_observer, - obj, - self.other_observer, - other_freq * 1e-9, - self.other_unit, + self.map_observer, obj, self.other_observer, other_freq * 1e-9, self.other_unit ), wrap=True, ) @@ -1339,9 +1154,7 @@ class overplot_radio(align_maps): (self.cr_map,) = self.ax_overplot.plot(*(self.map_wcs.celestial.wcs.crpix - (1.0, 1.0)), "r+") (self.cr_other,) = self.ax_overplot.plot( - *(self.other_wcs.celestial.wcs.crpix - (1.0, 1.0)), - "g+", - transform=self.ax_overplot.get_transform(self.other_wcs), + *(self.other_wcs.celestial.wcs.crpix - (1.0, 1.0)), "g+", transform=self.ax_overplot.get_transform(self.other_wcs) ) handles, labels = self.ax_overplot.get_legend_handles_labels() @@ -1351,12 +1164,7 @@ class overplot_radio(align_maps): labels.append("{0:s} contour".format(self.other_observer)) handles.append(Rectangle((0, 0), 1, 1, fill=False, lw=2, ec=other_cont.get_edgecolor()[0])) self.legend = self.ax_overplot.legend( - handles=handles, - labels=labels, - bbox_to_anchor=(0.0, 1.02, 1.0, 0.102), - loc="lower left", - mode="expand", - borderaxespad=0.0, + handles=handles, labels=labels, bbox_to_anchor=(0.0, 1.02, 1.0, 0.102), loc="lower left", mode="expand", borderaxespad=0.0 ) if savename is not None: @@ -1379,16 +1187,7 @@ class overplot_chandra(align_maps): Inherit from class align_maps in order to get the same WCS on both maps. """ - def overplot( - self, - levels=None, - P_cut=0.99, - SNRi_cut=1.0, - scale_vec=2, - zoom=1, - savename=None, - **kwargs, - ): + def overplot(self, levels=None, P_cut=0.99, SNRi_cut=1.0, scale_vec=2, zoom=1, savename=None, **kwargs): self.Stokes_UV = self.map self.wcs_UV = self.map_wcs # Get Data @@ -1402,10 +1201,7 @@ class overplot_chandra(align_maps): pang = self.Stokes_UV["POL_ANG"].data # Compute confidence level map QN, UN, QN_ERR, UN_ERR = np.full((4, stkI.shape[0], stkI.shape[1]), np.nan) - for nflux, sflux in zip( - [QN, UN, QN_ERR, UN_ERR], - [stkQ, stkU, np.sqrt(stk_cov[1, 1]), np.sqrt(stk_cov[2, 2])], - ): + for nflux, sflux in zip([QN, UN, QN_ERR, UN_ERR], [stkQ, stkU, np.sqrt(stk_cov[1, 1]), np.sqrt(stk_cov[2, 2])]): nflux[stkI > 0.0] = sflux[stkI > 0.0] / stkI[stkI > 0.0] confP = PCconf(QN, UN, QN_ERR, UN_ERR) @@ -1438,14 +1234,8 @@ class overplot_chandra(align_maps): self.fig_overplot.subplots_adjust(hspace=0, wspace=0, bottom=0.1, left=0.1, top=0.8, right=1) # Display UV intensity map with polarization vectors - vmin, vmax = ( - stkI[np.isfinite(stkI)].max() / 1e3 * self.map_convert, - stkI[np.isfinite(stkI)].max() * self.map_convert, - ) - for key, value in [ - ["cmap", [["cmap", "inferno"]]], - ["norm", [["norm", LogNorm(vmin, vmax)]]], - ]: + vmin, vmax = (stkI[np.isfinite(stkI)].max() / 1e3 * self.map_convert, stkI[np.isfinite(stkI)].max() * self.map_convert) + for key, value in [["cmap", [["cmap", "inferno"]]], ["norm", [["norm", LogNorm(vmin, vmax)]]]]: try: _ = kwargs[key] except KeyError: @@ -1482,12 +1272,7 @@ class overplot_chandra(align_maps): font_color = "black" self.im = self.ax_overplot.imshow(stkI * self.map_convert, aspect="equal", **kwargs) self.cbar = self.fig_overplot.colorbar( - self.im, - ax=self.ax_overplot, - aspect=50, - shrink=0.75, - pad=0.025, - label=r"$F_{{\lambda}}$ [{0:s}]".format(self.map_unit), + self.im, ax=self.ax_overplot, aspect=50, shrink=0.75, pad=0.025, label=r"$F_{{\lambda}}$ [{0:s}]".format(self.map_unit) ) # Display full size polarization vectors @@ -1498,10 +1283,7 @@ class overplot_chandra(align_maps): self.scale_vec = scale_vec step_vec = 1 self.X, self.Y = np.meshgrid(np.arange(stkI.shape[1]), np.arange(stkI.shape[0])) - self.U, self.V = ( - pol * np.cos(np.pi / 2.0 + pang * np.pi / 180.0), - pol * np.sin(np.pi / 2.0 + pang * np.pi / 180.0), - ) + self.U, self.V = (pol * np.cos(np.pi / 2.0 + pang * np.pi / 180.0), pol * np.sin(np.pi / 2.0 + pang * np.pi / 180.0)) self.Q = self.ax_overplot.quiver( self.X[::step_vec, ::step_vec], self.Y[::step_vec, ::step_vec], @@ -1529,18 +1311,14 @@ class overplot_chandra(align_maps): elif zoom != 1: levels *= other_data.max() / self.other_data.max() other_cont = self.ax_overplot.contour( - other_data * self.other_convert, - transform=self.ax_overplot.get_transform(other_wcs), - levels=levels, - colors="grey", + other_data * self.other_convert, transform=self.ax_overplot.get_transform(other_wcs), levels=levels, colors="grey" ) self.ax_overplot.clabel(other_cont, inline=True, fontsize=8) self.ax_overplot.set_xlabel(label="Right Ascension (J2000)") self.ax_overplot.set_ylabel(label="Declination (J2000)", labelpad=-1) self.fig_overplot.suptitle( - "{0:s} polarization map of {1:s} overplotted\nwith {2:s} contour in counts.".format(self.map_observer, obj, self.other_observer), - wrap=True, + "{0:s} polarization map of {1:s} overplotted\nwith {2:s} contour in counts.".format(self.map_observer, obj, self.other_observer), wrap=True ) # Display pixel scale and North direction @@ -1591,11 +1369,7 @@ class overplot_chandra(align_maps): self.ax_overplot.add_artist(pol_sc) (self.cr_map,) = self.ax_overplot.plot(*(self.map_wcs.celestial.wcs.crpix - (1.0, 1.0)), "r+") - (self.cr_other,) = self.ax_overplot.plot( - *(other_wcs.celestial.wcs.crpix - (1.0, 1.0)), - "g+", - transform=self.ax_overplot.get_transform(other_wcs), - ) + (self.cr_other,) = self.ax_overplot.plot(*(other_wcs.celestial.wcs.crpix - (1.0, 1.0)), "g+", transform=self.ax_overplot.get_transform(other_wcs)) handles, labels = self.ax_overplot.get_legend_handles_labels() handles[np.argmax([li == "{0:s} polarization map".format(self.map_observer) for li in labels])] = FancyArrowPatch( (0, 0), (0, 1), arrowstyle="-", fc="w", ec="k", lw=2 @@ -1603,12 +1377,7 @@ class overplot_chandra(align_maps): labels.append("{0:s} contour in counts".format(self.other_observer)) handles.append(Rectangle((0, 0), 1, 1, fill=False, lw=2, ec=other_cont.get_edgecolor()[0])) self.legend = self.ax_overplot.legend( - handles=handles, - labels=labels, - bbox_to_anchor=(0.0, 1.02, 1.0, 0.102), - loc="lower left", - mode="expand", - borderaxespad=0.0, + handles=handles, labels=labels, bbox_to_anchor=(0.0, 1.02, 1.0, 0.102), loc="lower left", mode="expand", borderaxespad=0.0 ) if savename is not None: @@ -1621,14 +1390,7 @@ class overplot_chandra(align_maps): def plot(self, levels=None, P_cut=0.99, SNRi_cut=1.0, zoom=1, savename=None, **kwargs) -> None: while not self.aligned: self.align() - self.overplot( - levels=levels, - P_cut=P_cut, - SNRi_cut=SNRi_cut, - zoom=zoom, - savename=savename, - **kwargs, - ) + self.overplot(levels=levels, P_cut=P_cut, SNRi_cut=SNRi_cut, zoom=zoom, savename=savename, **kwargs) plt.show(block=True) @@ -1638,17 +1400,7 @@ class overplot_pol(align_maps): Inherit from class align_maps in order to get the same WCS on both maps. """ - def overplot( - self, - levels="Default", - P_cut=0.99, - SNRi_cut=1.0, - step_vec=1, - scale_vec=2.0, - disptype="i", - savename=None, - **kwargs, - ): + def overplot(self, levels="Default", P_cut=0.99, SNRi_cut=1.0, step_vec=1, scale_vec=2.0, disptype="i", savename=None, **kwargs): self.Stokes_UV = self.map self.wcs_UV = self.map_wcs # Get Data @@ -1662,10 +1414,7 @@ class overplot_pol(align_maps): pang = self.Stokes_UV["POL_ANG"].data # Compute confidence level map QN, UN, QN_ERR, UN_ERR = np.full((4, stkI.shape[0], stkI.shape[1]), np.nan) - for nflux, sflux in zip( - [QN, UN, QN_ERR, UN_ERR], - [stkQ, stkU, np.sqrt(stk_cov[1, 1]), np.sqrt(stk_cov[2, 2])], - ): + for nflux, sflux in zip([QN, UN, QN_ERR, UN_ERR], [stkQ, stkU, np.sqrt(stk_cov[1, 1]), np.sqrt(stk_cov[2, 2])]): nflux[stkI > 0.0] = sflux[stkI > 0.0] / stkI[stkI > 0.0] confP = PCconf(QN, UN, QN_ERR, UN_ERR) @@ -1689,19 +1438,13 @@ class overplot_pol(align_maps): ratiox = max(int(stkI.shape[1] / stkI.shape[0]), 1) ratioy = max(int(stkI.shape[0] / stkI.shape[1]), 1) self.px_scale = self.wcs_UV.wcs.get_cdelt()[0] / self.other_wcs.wcs.get_cdelt()[0] - self.fig_overplot, self.ax_overplot = plt.subplots( - figsize=(10 * ratiox, 12 * ratioy), - subplot_kw=dict(projection=self.other_wcs), - ) + self.fig_overplot, self.ax_overplot = plt.subplots(figsize=(10 * ratiox, 12 * ratioy), subplot_kw=dict(projection=self.other_wcs)) self.fig_overplot.subplots_adjust(hspace=0, wspace=0, bottom=0.1, left=0.1, top=0.80, right=1.02) self.ax_overplot.set_xlabel(label="Right Ascension (J2000)") self.ax_overplot.set_ylabel(label="Declination (J2000)", labelpad=-1) # Display "other" intensity map - vmin, vmax = ( - other_data[other_data > 0.0].max() / 1e3 * self.other_convert, - other_data[other_data > 0.0].max() * self.other_convert, - ) + vmin, vmax = (other_data[other_data > 0.0].max() / 1e3 * self.other_convert, other_data[other_data > 0.0].max() * self.other_convert) for key, value in [ ["cmap", [["cmap", "inferno"]]], ["norm", [["vmin", vmin], ["vmax", vmax]]], @@ -1744,11 +1487,7 @@ class overplot_pol(align_maps): font_color = "black" if "norm" in kwargs.keys(): self.im = self.ax_overplot.imshow( - other_data * self.other_convert, - alpha=1.0, - label="{0:s} observation".format(self.other_observer), - cmap=kwargs["cmap"], - norm=kwargs["norm"], + other_data * self.other_convert, alpha=1.0, label="{0:s} observation".format(self.other_observer), cmap=kwargs["cmap"], norm=kwargs["norm"] ) else: self.im = self.ax_overplot.imshow( @@ -1760,12 +1499,7 @@ class overplot_pol(align_maps): vmax=kwargs["vmax"], ) self.cbar = self.fig_overplot.colorbar( - self.im, - ax=self.ax_overplot, - aspect=80, - shrink=0.75, - pad=0.025, - label=r"$F_{{\lambda}}$ [{0:s}]".format(self.other_unit), + self.im, ax=self.ax_overplot, aspect=80, shrink=0.75, pad=0.025, label=r"$F_{{\lambda}}$ [{0:s}]".format(self.other_unit) ) # Display full size polarization vectors @@ -1778,10 +1512,7 @@ class overplot_pol(align_maps): else: self.scale_vec = scale_vec * self.px_scale self.X, self.Y = np.meshgrid(np.arange(stkI.shape[1]), np.arange(stkI.shape[0])) - self.U, self.V = ( - pol * np.cos(np.pi / 2.0 + pang * np.pi / 180.0), - pol * np.sin(np.pi / 2.0 + pang * np.pi / 180.0), - ) + self.U, self.V = (pol * np.cos(np.pi / 2.0 + pang * np.pi / 180.0), pol * np.sin(np.pi / 2.0 + pang * np.pi / 180.0)) self.Q = self.ax_overplot.quiver( self.X[::step_vec, ::step_vec], self.Y[::step_vec, ::step_vec], @@ -1811,44 +1542,26 @@ class overplot_pol(align_maps): if levels == "Default": levels = np.array([2.0, 5.0, 10.0, 20.0, 90.0]) / 100.0 * np.max(pol[stkI > 0.0]) cont_stk = self.ax_overplot.contour( - pol * 100.0, - levels=levels * 100.0, - colors="grey", - alpha=0.75, - transform=self.ax_overplot.get_transform(self.wcs_UV), + pol * 100.0, levels=levels * 100.0, colors="grey", alpha=0.75, transform=self.ax_overplot.get_transform(self.wcs_UV) ) elif disptype.lower() == "pf": disptypestr = "polarized flux" if levels == "Default": levels = np.array([2.0, 5.0, 10.0, 20.0, 90.0]) / 100.0 * np.max(stkI[stkI > 0.0] * pol[stkI > 0.0]) * self.map_convert cont_stk = self.ax_overplot.contour( - stkI * pol * self.map_convert, - levels=levels, - colors="grey", - alpha=0.75, - transform=self.ax_overplot.get_transform(self.wcs_UV), + stkI * pol * self.map_convert, levels=levels, colors="grey", alpha=0.75, transform=self.ax_overplot.get_transform(self.wcs_UV) ) elif disptype.lower() == "snri": disptypestr = "Stokes I signal-to-noise" if levels == "Default": levels = np.array([2.0, 5.0, 10.0, 20.0, 90.0]) / 100.0 * np.max(SNRi[stk_cov[0, 0] > 0.0]) - cont_stk = self.ax_overplot.contour( - SNRi, - levels=levels, - colors="grey", - alpha=0.75, - transform=self.ax_overplot.get_transform(self.wcs_UV), - ) + cont_stk = self.ax_overplot.contour(SNRi, levels=levels, colors="grey", alpha=0.75, transform=self.ax_overplot.get_transform(self.wcs_UV)) else: # default to intensity contours disptypestr = "Stokes I" if levels == "Default": levels = np.array([2.0, 5.0, 10.0, 20.0, 90.0]) / 100.0 * np.max(stkI[stkI > 0.0]) * self.map_convert cont_stk = self.ax_overplot.contour( - stkI * self.map_convert, - levels=levels, - colors="grey", - alpha=0.75, - transform=self.ax_overplot.get_transform(self.wcs_UV), + stkI * self.map_convert, levels=levels, colors="grey", alpha=0.75, transform=self.ax_overplot.get_transform(self.wcs_UV) ) # self.ax_overplot.clabel(cont_stk, inline=False, colors="k", fontsize=7) @@ -1902,11 +1615,7 @@ class overplot_pol(align_maps): ) self.ax_overplot.add_artist(pol_sc) - (self.cr_map,) = self.ax_overplot.plot( - *(self.map_wcs.celestial.wcs.crpix - (1.0, 1.0)), - "r+", - transform=self.ax_overplot.get_transform(self.wcs_UV), - ) + (self.cr_map,) = self.ax_overplot.plot(*(self.map_wcs.celestial.wcs.crpix - (1.0, 1.0)), "r+", transform=self.ax_overplot.get_transform(self.wcs_UV)) (self.cr_other,) = self.ax_overplot.plot(*(self.other_wcs.celestial.wcs.crpix - (1.0, 1.0)), "g+") if "PHOTPLAM" in list(self.other_header.keys()): @@ -1926,12 +1635,7 @@ class overplot_pol(align_maps): handles.append(Rectangle((0, 0), 1, 1, fill=False, ec=cont_stk.get_edgecolor()[0])) disptypestr += " contours" self.legend = self.ax_overplot.legend( - handles=handles, - labels=labels, - bbox_to_anchor=(0.0, 1.02, 1.0, 0.102), - loc="lower left", - mode="expand", - borderaxespad=0.0, + handles=handles, labels=labels, bbox_to_anchor=(0.0, 1.02, 1.0, 0.102), loc="lower left", mode="expand", borderaxespad=0.0 ) self.fig_overplot.suptitle( @@ -1946,29 +1650,10 @@ class overplot_pol(align_maps): self.fig_overplot.canvas.draw() - def plot( - self, - levels=None, - P_cut=0.99, - SNRi_cut=1.0, - step_vec=1, - scale_vec=2.0, - disptype="i", - savename=None, - **kwargs, - ) -> None: + def plot(self, levels=None, P_cut=0.99, SNRi_cut=1.0, step_vec=1, scale_vec=2.0, disptype="i", savename=None, **kwargs) -> None: while not self.aligned: self.align() - self.overplot( - levels=levels, - P_cut=P_cut, - SNRi_cut=SNRi_cut, - step_vec=step_vec, - scale_vec=scale_vec, - disptype=disptype, - savename=savename, - **kwargs, - ) + self.overplot(levels=levels, P_cut=P_cut, SNRi_cut=SNRi_cut, step_vec=step_vec, scale_vec=scale_vec, disptype=disptype, savename=savename, **kwargs) plt.show(block=True) def add_vector(self, position="center", pol_deg=1.0, pol_ang=0.0, **kwargs): @@ -1977,10 +1662,7 @@ class overplot_pol(align_maps): if isinstance(position, SkyCoord): position = self.other_wcs.world_to_pixel(position) - u, v = ( - pol_deg * np.cos(np.radians(pol_ang) + np.pi / 2.0), - pol_deg * np.sin(np.radians(pol_ang) + np.pi / 2.0), - ) + u, v = (pol_deg * np.cos(np.radians(pol_ang) + np.pi / 2.0), pol_deg * np.sin(np.radians(pol_ang) + np.pi / 2.0)) for key, value in [ ["scale", [["scale", 1.0 / self.scale_vec]]], ["width", [["width", 0.5 * self.px_scale]]], @@ -1993,34 +1675,15 @@ class overplot_pol(align_maps): except KeyError: for key_i, val_i in value: kwargs[key_i] = val_i - handles, labels = ( - self.legend.legend_handles, - [l.get_text() for l in self.legend.texts], - ) + handles, labels = (self.legend.legend_handles, [l.get_text() for l in self.legend.texts]) new_vec = self.ax_overplot.quiver( - *position, - u, - v, - units="xy", - angles="uv", - scale_units="xy", - pivot="mid", - headwidth=0.0, - headlength=0.0, - headaxislength=0.0, - **kwargs, + *position, u, v, units="xy", angles="uv", scale_units="xy", pivot="mid", headwidth=0.0, headlength=0.0, headaxislength=0.0, **kwargs ) handles.append(FancyArrowPatch((0, 0), (0, 1), arrowstyle="-", fc=new_vec.get_fc(), ec=new_vec.get_ec())) labels.append(new_vec.get_label()) self.legend.remove() self.legend = self.ax_overplot.legend( - title=self.legend_title, - handles=handles, - labels=labels, - bbox_to_anchor=(0.0, 1.02, 1.0, 0.102), - loc="lower left", - mode="expand", - borderaxespad=0.0, + title=self.legend_title, handles=handles, labels=labels, bbox_to_anchor=(0.0, 1.02, 1.0, 0.102), loc="lower left", mode="expand", borderaxespad=0.0 ) self.fig_overplot.canvas.draw() return new_vec @@ -2039,17 +1702,7 @@ class align_pol(object): self.kwargs = kwargs - def single_plot( - self, - curr_map, - wcs, - v_lim=None, - ax_lim=None, - P_cut=3.0, - SNRi_cut=3.0, - savename=None, - **kwargs, - ): + def single_plot(self, curr_map, wcs, v_lim=None, ax_lim=None, P_cut=3.0, SNRi_cut=3.0, savename=None, **kwargs): # Get data stkI = curr_map["I_STOKES"].data stk_cov = curr_map["IQU_COV_MATRIX"].data @@ -2098,10 +1751,7 @@ class align_pol(object): else: vmin, vmax = v_lim * convert_flux - for key, value in [ - ["cmap", [["cmap", "inferno"]]], - ["norm", [["vmin", vmin], ["vmax", vmax]]], - ]: + for key, value in [["cmap", [["cmap", "inferno"]]], ["norm", [["vmin", vmin], ["vmax", vmax]]]]: try: test = kwargs[key] if isinstance(test, LogNorm): @@ -2111,28 +1761,10 @@ class align_pol(object): kwargs[key_i] = val_i im = ax.imshow(stkI * convert_flux, aspect="equal", **kwargs) - fig.colorbar( - im, - ax=ax, - aspect=50, - shrink=0.75, - pad=0.025, - label=r"$F_{\lambda}$ [$ergs \cdot cm^{-2} \cdot s^{-1} \cdot \AA^{-1}$]", - ) + fig.colorbar(im, ax=ax, aspect=50, shrink=0.75, pad=0.025, label=r"$F_{\lambda}$ [$ergs \cdot cm^{-2} \cdot s^{-1} \cdot \AA^{-1}$]") px_size = wcs.wcs.get_cdelt()[0] * 3600.0 - px_sc = AnchoredSizeBar( - ax.transData, - 1.0 / px_size, - "1 arcsec", - 3, - pad=0.5, - sep=5, - borderpad=0.5, - frameon=False, - size_vertical=0.005, - color="w", - ) + px_sc = AnchoredSizeBar(ax.transData, 1.0 / px_size, "1 arcsec", 3, pad=0.5, sep=5, borderpad=0.5, frameon=False, size_vertical=0.005, color="w") ax.add_artist(px_sc) north_dir = AnchoredDirectionArrows( @@ -2157,10 +1789,7 @@ class align_pol(object): step_vec = 1 X, Y = np.meshgrid(np.arange(stkI.shape[1]), np.arange(stkI.shape[0])) - U, V = ( - pol * np.cos(np.pi / 2.0 + pang * np.pi / 180.0), - pol * np.sin(np.pi / 2.0 + pang * np.pi / 180.0), - ) + U, V = (pol * np.cos(np.pi / 2.0 + pang * np.pi / 180.0), pol * np.sin(np.pi / 2.0 + pang * np.pi / 180.0)) ax.quiver( X[::step_vec, ::step_vec], Y[::step_vec, ::step_vec], @@ -2178,18 +1807,7 @@ class align_pol(object): linewidth=0.75, color="w", ) - pol_sc = AnchoredSizeBar( - ax.transData, - 2.0, - r"$P$= 100 %", - 4, - pad=0.5, - sep=5, - borderpad=0.5, - frameon=False, - size_vertical=0.005, - color="w", - ) + pol_sc = AnchoredSizeBar(ax.transData, 2.0, r"$P$= 100 %", 4, pad=0.5, sep=5, borderpad=0.5, frameon=False, size_vertical=0.005, color="w") ax.add_artist(pol_sc) if "PHOTPLAM" in list(curr_map[0].header.keys()): @@ -2224,17 +1842,7 @@ class align_pol(object): np.min( [ np.min( - curr_map[0].data[ - curr_map[0].data - > SNRi_cut - * np.max( - [ - eps * np.ones(curr_map[0].data.shape), - np.sqrt(curr_map[3].data[0, 0]), - ], - axis=0, - ) - ] + curr_map[0].data[curr_map[0].data > SNRi_cut * np.max([eps * np.ones(curr_map[0].data.shape), np.sqrt(curr_map[3].data[0, 0])], axis=0)] ) for curr_map in self.other_maps ] @@ -2243,19 +1851,7 @@ class align_pol(object): ) vmax = np.max( [ - np.max( - curr_map[0].data[ - curr_map[0].data - > SNRi_cut - * np.max( - [ - eps * np.ones(curr_map[0].data.shape), - np.sqrt(curr_map[3].data[0, 0]), - ], - axis=0, - ) - ] - ) + np.max(curr_map[0].data[curr_map[0].data > SNRi_cut * np.max([eps * np.ones(curr_map[0].data.shape), np.sqrt(curr_map[3].data[0, 0])], axis=0)]) for curr_map in self.other_maps ] ) @@ -2265,15 +1861,7 @@ class align_pol(object): vmin, np.min( self.ref_map[0].data[ - self.ref_map[0].data - > SNRi_cut - * np.max( - [ - eps * np.ones(self.ref_map[0].data.shape), - np.sqrt(self.ref_map[3].data[0, 0]), - ], - axis=0, - ) + self.ref_map[0].data > SNRi_cut * np.max([eps * np.ones(self.ref_map[0].data.shape), np.sqrt(self.ref_map[3].data[0, 0])], axis=0) ] ), ] @@ -2285,43 +1873,20 @@ class align_pol(object): vmax, np.max( self.ref_map[0].data[ - self.ref_map[0].data - > SNRi_cut - * np.max( - [ - eps * np.ones(self.ref_map[0].data.shape), - np.sqrt(self.ref_map[3].data[0, 0]), - ], - axis=0, - ) + self.ref_map[0].data > SNRi_cut * np.max([eps * np.ones(self.ref_map[0].data.shape), np.sqrt(self.ref_map[3].data[0, 0])], axis=0) ] ), ] ) v_lim = np.array([vmin, vmax]) - fig, ax = self.single_plot( - self.ref_map, - self.wcs, - v_lim=v_lim, - P_cut=P_cut, - SNRi_cut=SNRi_cut, - savename=savename + "_0", - **kwargs, - ) + fig, ax = self.single_plot(self.ref_map, self.wcs, v_lim=v_lim, P_cut=P_cut, SNRi_cut=SNRi_cut, savename=savename + "_0", **kwargs) x_lim, y_lim = ax.get_xlim(), ax.get_ylim() ax_lim = np.array([self.wcs.pixel_to_world(x_lim[i], y_lim[i]) for i in range(len(x_lim))]) for i, curr_map in enumerate(self.other_maps): self.single_plot( - curr_map, - self.wcs_other[i], - v_lim=v_lim, - ax_lim=ax_lim, - P_cut=P_cut, - SNRi_cut=SNRi_cut, - savename=savename + "_" + str(i + 1), - **kwargs, + curr_map, self.wcs_other[i], v_lim=v_lim, ax_lim=ax_lim, P_cut=P_cut, SNRi_cut=SNRi_cut, savename=savename + "_" + str(i + 1), **kwargs ) @@ -2388,10 +1953,7 @@ class crop_map(object): else: kwargs = {**self.kwargs, **kwargs} - vmin, vmax = ( - np.min(data[data > 0.0] * convert_flux), - np.max(data[data > 0.0] * convert_flux), - ) + vmin, vmax = (np.min(data[data > 0.0] * convert_flux), np.max(data[data > 0.0] * convert_flux)) for key, value in [ ["cmap", [["cmap", "inferno"]]], ["origin", [["origin", "lower"]]], @@ -2604,15 +2166,9 @@ class crop_Stokes(crop_map): if dataset.header["FILENAME"][-4:] != "crop": dataset.header["FILENAME"] += "_crop" dataset.header["P_int"] = (P_diluted, "Integrated polarization degree") - dataset.header["sP_int"] = ( - np.ceil(P_diluted_err * 1000.0) / 1000.0, - "Integrated polarization degree error", - ) + dataset.header["sP_int"] = (np.ceil(P_diluted_err * 1000.0) / 1000.0, "Integrated polarization degree error") dataset.header["PA_int"] = (PA_diluted, "Integrated polarization angle") - dataset.header["sPA_int"] = ( - np.ceil(PA_diluted_err * 10.0) / 10.0, - "Integrated polarization angle error", - ) + dataset.header["sPA_int"] = (np.ceil(PA_diluted_err * 10.0) / 10.0, "Integrated polarization angle error") self.fig.canvas.draw_idle() @property @@ -2641,14 +2197,7 @@ class image_lasso_selector(object): self.ax = ax self.mask_alpha = 0.1 self.embedded = True - self.displayed = self.ax.imshow( - self.img, - vmin=self.vmin, - vmax=self.vmax, - aspect="equal", - cmap="inferno", - alpha=self.mask_alpha, - ) + self.displayed = self.ax.imshow(self.img, vmin=self.vmin, vmax=self.vmax, aspect="equal", cmap="inferno", alpha=self.mask_alpha) plt.ion() lineprops = {"color": "grey", "linewidth": 1, "alpha": 0.8} @@ -2677,14 +2226,7 @@ class image_lasso_selector(object): def update_mask(self): self.displayed.remove() - self.displayed = self.ax.imshow( - self.img, - vmin=self.vmin, - vmax=self.vmax, - aspect="equal", - cmap="inferno", - alpha=self.mask_alpha, - ) + self.displayed = self.ax.imshow(self.img, vmin=self.vmin, vmax=self.vmax, aspect="equal", cmap="inferno", alpha=self.mask_alpha) array = self.displayed.get_array().data self.mask = np.zeros(self.img.shape[:2], dtype=bool) @@ -2700,16 +2242,7 @@ class image_lasso_selector(object): class slit(object): - def __init__( - self, - img, - cdelt=np.array([1.0, 1.0]), - width=1.0, - height=2.0, - angle=0.0, - fig=None, - ax=None, - ): + def __init__(self, img, cdelt=np.array([1.0, 1.0]), width=1.0, height=2.0, angle=0.0, fig=None, ax=None): """ img must have shape (X, Y) """ @@ -2730,14 +2263,7 @@ class slit(object): self.mask_alpha = 0.1 self.embedded = True - self.displayed = self.ax.imshow( - self.img, - vmin=self.vmin, - vmax=self.vmax, - aspect="equal", - cmap="inferno", - alpha=self.mask_alpha, - ) + self.displayed = self.ax.imshow(self.img, vmin=self.vmin, vmax=self.vmax, aspect="equal", cmap="inferno", alpha=self.mask_alpha) plt.ion() xx, yy = np.indices(self.img.shape) @@ -2751,15 +2277,7 @@ class slit(object): self.angle = angle self.rect_center = (self.x0, self.y0) - np.dot(rot2D(self.angle), (self.width / 2, self.height / 2)) - self.rect = Rectangle( - self.rect_center, - self.width, - self.height, - angle=self.angle, - alpha=0.8, - ec="grey", - fc="none", - ) + self.rect = Rectangle(self.rect_center, self.width, self.height, angle=self.angle, alpha=0.8, ec="grey", fc="none") self.ax.add_patch(self.rect) self.fig.canvas.mpl_connect("button_press_event", self.on_press) @@ -2819,14 +2337,7 @@ class slit(object): self.displayed.remove() except ValueError: return - self.displayed = self.ax.imshow( - self.img, - vmin=self.vmin, - vmax=self.vmax, - aspect="equal", - cmap="inferno", - alpha=self.mask_alpha, - ) + self.displayed = self.ax.imshow(self.img, vmin=self.vmin, vmax=self.vmax, aspect="equal", cmap="inferno", alpha=self.mask_alpha) array = self.displayed.get_array().data self.mask = np.zeros(array.shape, dtype=bool) @@ -2864,14 +2375,7 @@ class aperture(object): self.mask_alpha = 0.1 self.embedded = True - self.displayed = self.ax.imshow( - self.img, - vmin=self.vmin, - vmax=self.vmax, - aspect="equal", - cmap="inferno", - alpha=self.mask_alpha, - ) + self.displayed = self.ax.imshow(self.img, vmin=self.vmin, vmax=self.vmax, aspect="equal", cmap="inferno", alpha=self.mask_alpha) plt.ion() xx, yy = np.indices(self.img.shape) @@ -2932,14 +2436,7 @@ class aperture(object): self.displayed.remove() except ValueError: return - self.displayed = self.ax.imshow( - self.img, - vmin=self.vmin, - vmax=self.vmax, - aspect="equal", - cmap="inferno", - alpha=self.mask_alpha, - ) + self.displayed = self.ax.imshow(self.img, vmin=self.vmin, vmax=self.vmax, aspect="equal", cmap="inferno", alpha=self.mask_alpha) array = self.displayed.get_array().data yy, xx = np.indices(self.img.shape[:2]) @@ -2960,17 +2457,7 @@ class pol_map(object): Class to interactively study polarization maps. """ - def __init__( - self, - Stokes, - P_cut=0.99, - SNRi_cut=1.0, - step_vec=1, - scale_vec=3.0, - flux_lim=None, - selection=None, - pa_err=False, - ): + def __init__(self, Stokes, P_cut=0.99, SNRi_cut=1.0, step_vec=1, scale_vec=3.0, flux_lim=None, selection=None, pa_err=False): if isinstance(Stokes, str): Stokes = fits.open(Stokes) self.Stokes = deepcopy(Stokes) @@ -3013,22 +2500,8 @@ class pol_map(object): ax_snr_conf = self.fig.add_axes([0.115, 0.020, 0.05, 0.02]) SNRi_max = np.max(self.I[self.IQU_cov[0, 0] > 0.0] / np.sqrt(self.IQU_cov[0, 0][self.IQU_cov[0, 0] > 0.0])) SNRp_max = np.max(self.P[self.s_P > 0.0] / self.s_P[self.s_P > 0.0]) - s_I_cut = Slider( - ax_I_cut, - r"$SNR^{I}_{cut}$", - 1.0, - int(SNRi_max * 0.95), - valstep=1, - valinit=self.SNRi_cut, - ) - self.s_P_cut = Slider( - self.ax_P_cut, - r"$Conf^{P}_{cut}$", - 0.50, - 1.0, - valstep=[0.500, 0.900, 0.990, 0.999], - valinit=self.P_cut if P_cut <= 1.0 else 0.99, - ) + s_I_cut = Slider(ax_I_cut, r"$SNR^{I}_{cut}$", 1.0, int(SNRi_max * 0.95), valstep=1, valinit=self.SNRi_cut) + self.s_P_cut = Slider(self.ax_P_cut, r"$Conf^{P}_{cut}$", 0.50, 1.0, valstep=[0.500, 0.900, 0.990, 0.999], valinit=self.P_cut if P_cut <= 1.0 else 0.99) s_vec_sc = Slider(ax_vec_sc, r"Vec scale", 0.0, 10.0, valstep=1, valinit=self.scale_vec) b_snr_reset = Button(ax_snr_reset, "Reset") b_snr_reset.label.set_fontsize(8) @@ -3072,12 +2545,7 @@ class pol_map(object): self.snr_conf = 1 b_snr_conf.label.set_text("SNR") self.s_P_cut = Slider( - self.ax_P_cut, - r"$Conf^{P}_{cut}$", - 0.50, - 1.0, - valstep=[0.500, 0.900, 0.990, 0.999], - valinit=self.P_cut if P_cut <= 1.0 else 0.99, + self.ax_P_cut, r"$Conf^{P}_{cut}$", 0.50, 1.0, valstep=[0.500, 0.900, 0.990, 0.999], valinit=self.P_cut if P_cut <= 1.0 else 0.99 ) self.s_P_cut.on_changed(update_snrp) update_snrp(self.s_P_cut.val) @@ -3141,14 +2609,7 @@ class pol_map(object): b_aper.label.set_fontsize(8) b_aper_reset = Button(ax_aper_reset, "Reset") b_aper_reset.label.set_fontsize(8) - s_aper_radius = Slider( - ax_aper_radius, - r"$R_{aper}$", - np.ceil(self.wcs.wcs.cdelt.max() / 1.33 * 3.6e5) / 1e2, - 3.5, - valstep=1e-2, - valinit=1.0, - ) + s_aper_radius = Slider(ax_aper_radius, r"$R_{aper}$", np.ceil(self.wcs.wcs.cdelt.max() / 1.33 * 3.6e5) / 1e2, 3.5, valstep=1e-2, valinit=1.0) def select_aperture(event): if self.data is None: @@ -3165,13 +2626,7 @@ class pol_map(object): else: self.selected = True self.region = None - self.select_instance = aperture( - self.data, - fig=self.fig, - ax=self.ax, - cdelt=self.wcs.wcs.cdelt, - radius=s_aper_radius.val, - ) + self.select_instance = aperture(self.data, fig=self.fig, ax=self.ax, cdelt=self.wcs.wcs.cdelt, radius=s_aper_radius.val) self.select_instance.circ.set_visible(True) self.fig.canvas.draw_idle() @@ -3182,22 +2637,10 @@ class pol_map(object): self.select_instance.update_radius(val) else: self.selected = True - self.select_instance = aperture( - self.data, - fig=self.fig, - ax=self.ax, - cdelt=self.wcs.wcs.cdelt, - radius=val, - ) + self.select_instance = aperture(self.data, fig=self.fig, ax=self.ax, cdelt=self.wcs.wcs.cdelt, radius=val) else: self.selected = True - self.select_instance = aperture( - self.data, - fig=self.fig, - ax=self.ax, - cdelt=self.wcs.wcs.cdelt, - radius=val, - ) + self.select_instance = aperture(self.data, fig=self.fig, ax=self.ax, cdelt=self.wcs.wcs.cdelt, radius=val) self.fig.canvas.draw_idle() def reset_aperture(event): @@ -3221,22 +2664,8 @@ class pol_map(object): b_slit.label.set_fontsize(8) b_slit_reset = Button(ax_slit_reset, "Reset") b_slit_reset.label.set_fontsize(8) - s_slit_width = Slider( - ax_slit_width, - r"$W_{slit}$", - np.ceil(self.wcs.wcs.cdelt.max() / 1.33 * 3.6e5) / 1e2, - 7.0, - valstep=1e-2, - valinit=1.0, - ) - s_slit_height = Slider( - ax_slit_height, - r"$H_{slit}$", - np.ceil(self.wcs.wcs.cdelt.max() / 1.33 * 3.6e5) / 1e2, - 7.0, - valstep=1e-2, - valinit=1.0, - ) + s_slit_width = Slider(ax_slit_width, r"$W_{slit}$", np.ceil(self.wcs.wcs.cdelt.max() / 1.33 * 3.6e5) / 1e2, 7.0, valstep=1e-2, valinit=1.0) + s_slit_height = Slider(ax_slit_height, r"$H_{slit}$", np.ceil(self.wcs.wcs.cdelt.max() / 1.33 * 3.6e5) / 1e2, 7.0, valstep=1e-2, valinit=1.0) s_slit_angle = Slider(ax_slit_angle, r"$\theta_{slit}$", 0.0, 90.0, valstep=1.0, valinit=0.0) def select_slit(event): @@ -3255,13 +2684,7 @@ class pol_map(object): self.selected = True self.region = None self.select_instance = slit( - self.data, - fig=self.fig, - ax=self.ax, - cdelt=self.wcs.wcs.cdelt, - width=s_slit_width.val, - height=s_slit_height.val, - angle=s_slit_angle.val, + self.data, fig=self.fig, ax=self.ax, cdelt=self.wcs.wcs.cdelt, width=s_slit_width.val, height=s_slit_height.val, angle=s_slit_angle.val ) self.select_instance.rect.set_visible(True) @@ -3274,24 +2697,12 @@ class pol_map(object): else: self.selected = True self.select_instance = slit( - self.data, - fig=self.fig, - ax=self.ax, - cdelt=self.wcs.wcs.cdelt, - width=val, - height=s_slit_height.val, - angle=s_slit_angle.val, + self.data, fig=self.fig, ax=self.ax, cdelt=self.wcs.wcs.cdelt, width=val, height=s_slit_height.val, angle=s_slit_angle.val ) else: self.selected = True self.select_instance = slit( - self.data, - fig=self.fig, - ax=self.ax, - cdelt=self.wcs.wcs.cdelt, - width=val, - height=s_slit_height.val, - angle=s_slit_angle.val, + self.data, fig=self.fig, ax=self.ax, cdelt=self.wcs.wcs.cdelt, width=val, height=s_slit_height.val, angle=s_slit_angle.val ) self.fig.canvas.draw_idle() @@ -3302,24 +2713,12 @@ class pol_map(object): else: self.selected = True self.select_instance = slit( - self.data, - fig=self.fig, - ax=self.ax, - cdelt=self.wcs.wcs.cdelt, - width=s_slit_width.val, - height=val, - angle=s_slit_angle.val, + self.data, fig=self.fig, ax=self.ax, cdelt=self.wcs.wcs.cdelt, width=s_slit_width.val, height=val, angle=s_slit_angle.val ) else: self.selected = True self.select_instance = slit( - self.data, - fig=self.fig, - ax=self.ax, - cdelt=self.wcs.wcs.cdelt, - width=s_slit_width.val, - height=val, - angle=s_slit_angle.val, + self.data, fig=self.fig, ax=self.ax, cdelt=self.wcs.wcs.cdelt, width=s_slit_width.val, height=val, angle=s_slit_angle.val ) self.fig.canvas.draw_idle() @@ -3330,24 +2729,12 @@ class pol_map(object): else: self.selected = True self.select_instance = slit( - self.data, - fig=self.fig, - ax=self.ax, - cdelt=self.wcs.wcs.cdelt, - width=s_slit_width.val, - height=s_slit_height.val, - angle=val, + self.data, fig=self.fig, ax=self.ax, cdelt=self.wcs.wcs.cdelt, width=s_slit_width.val, height=s_slit_height.val, angle=val ) else: self.selected = True self.select_instance = slit( - self.data, - fig=self.fig, - ax=self.ax, - cdelt=self.wcs.wcs.cdelt, - width=s_slit_width.val, - height=s_slit_height.val, - angle=val, + self.data, fig=self.fig, ax=self.ax, cdelt=self.wcs.wcs.cdelt, width=s_slit_width.val, height=s_slit_height.val, angle=val ) self.fig.canvas.draw_idle() @@ -3432,11 +2819,7 @@ class pol_map(object): def submit_save(expression): ax_text_save.set(visible=False) if expression != "": - save_fig, save_ax = plt.subplots( - figsize=(12, 10), - layout="constrained", - subplot_kw=dict(projection=self.wcs), - ) + save_fig, save_ax = plt.subplots(figsize=(12, 10), layout="constrained", subplot_kw=dict(projection=self.wcs)) self.ax_cosmetics(ax=save_ax) self.display(fig=save_fig, ax=save_ax) self.pol_vector(fig=save_fig, ax=save_ax) @@ -3474,10 +2857,7 @@ class pol_map(object): center = (shape / 2).astype(int) cdelt_arcsec = self.wcs.wcs.cdelt * 3600 xx, yy = np.indices(shape) - x, y = ( - (xx - center[0]) * cdelt_arcsec[0], - (yy - center[1]) * cdelt_arcsec[1], - ) + x, y = ((xx - center[0]) * cdelt_arcsec[0], (yy - center[1]) * cdelt_arcsec[1]) P, PA = np.zeros(shape), np.zeros(shape) P[self.cut] = self.P[self.cut] @@ -3486,15 +2866,7 @@ class pol_map(object): for i in range(shape[0]): for j in range(shape[1]): dump_list.append( - [ - x[i, j], - y[i, j], - self.I[i, j] * self.map_convert, - self.Q[i, j] * self.map_convert, - self.U[i, j] * self.map_convert, - P[i, j], - PA[i, j], - ] + [x[i, j], y[i, j], self.I[i, j] * self.map_convert, self.Q[i, j] * self.map_convert, self.U[i, j] * self.map_convert, P[i, j], PA[i, j]] ) self.data_dump = np.array(dump_list) @@ -3659,10 +3031,7 @@ class pol_map(object): @property def cut(self): s_I = np.sqrt(self.IQU_cov[0, 0]) - SNRp_mask, SNRi_mask = ( - np.zeros(self.P.shape).astype(bool), - np.zeros(self.I.shape).astype(bool), - ) + SNRp_mask, SNRi_mask = (np.zeros(self.P.shape).astype(bool), np.zeros(self.I.shape).astype(bool)) SNRi_mask[s_I > 0.0] = self.I[s_I > 0.0] / s_I[s_I > 0.0] > self.SNRi if self.SNRp >= 1.0: SNRp_mask[self.s_P > 0.0] = self.P[self.s_P > 0.0] / self.s_P[self.s_P > 0.0] > self.SNRp @@ -3707,17 +3076,7 @@ class pol_map(object): else: scale_vec = 2.0 self.pol_sc = AnchoredSizeBar( - ax.transData, - scale_vec, - r"$P$= 100%", - 4, - pad=0.5, - sep=5, - borderpad=0.5, - frameon=False, - size_vertical=0.005, - color="white", - fontproperties=fontprops, + ax.transData, scale_vec, r"$P$= 100%", 4, pad=0.5, sep=5, borderpad=0.5, frameon=False, size_vertical=0.005, color="white", fontproperties=fontprops ) ax.add_artist(self.pol_sc) if hasattr(self, "north_dir"): @@ -3749,10 +3108,7 @@ class pol_map(object): if self.display_selection is None or self.display_selection.lower() in ["total_flux"]: self.data = self.Flux # self.I * self.map_convert if flux_lim is None: - vmin, vmax = ( - 1.0 / 2.0 * np.median(self.I[self.I > 0.0] * self.map_convert), - np.max(self.I[self.I > 0.0] * self.map_convert), - ) + vmin, vmax = (1.0 / 2.0 * np.median(self.I[self.I > 0.0] * self.map_convert), np.max(self.I[self.I > 0.0] * self.map_convert)) else: vmin, vmax = flux_lim kwargs["norm"] = LogNorm(vmin, vmax) @@ -3764,10 +3120,7 @@ class pol_map(object): elif self.display_selection.lower() in ["pol_flux"]: self.data = self.I * self.map_convert * self.P if flux_lim is None: - vmin, vmax = ( - 1.0 / 2.0 * np.median(self.I[self.I > 0.0] * self.map_convert), - np.max(self.I[self.I > 0.0] * self.map_convert), - ) + vmin, vmax = (1.0 / 2.0 * np.median(self.I[self.I > 0.0] * self.map_convert), np.max(self.I[self.I > 0.0] * self.map_convert)) else: vmin, vmax = flux_lim kwargs["norm"] = LogNorm(vmin, vmax) @@ -3828,10 +3181,7 @@ class pol_map(object): scale_vec = 2.0 P_cut[self.cut] = 1.0 / scale_vec X, Y = np.meshgrid(np.arange(self.I.shape[1]), np.arange(self.I.shape[0])) - XY_U, XY_V = ( - P_cut * np.cos(np.pi / 2.0 + self.PA * np.pi / 180.0), - P_cut * np.sin(np.pi / 2.0 + self.PA * np.pi / 180.0), - ) + XY_U, XY_V = (P_cut * np.cos(np.pi / 2.0 + self.PA * np.pi / 180.0), P_cut * np.sin(np.pi / 2.0 + self.PA * np.pi / 180.0)) if fig is None: fig = self.fig @@ -4062,12 +3412,7 @@ class pol_map(object): IU_reg_err = np.sqrt(np.sum(s_IU[self.region] ** 2)) QU_reg_err = np.sqrt(np.sum(s_QU[self.region] ** 2)) - conf = PCconf( - QN=Q_reg / I_reg, - QN_ERR=Q_reg_err / I_reg, - UN=U_reg / I_reg, - UN_ERR=U_reg_err / I_reg, - ) + conf = PCconf(QN=Q_reg / I_reg, QN_ERR=Q_reg_err / I_reg, UN=U_reg / I_reg, UN_ERR=U_reg_err / I_reg) if 1.0 - conf > 1e-3: str_conf = "\n" + r"Confidence: {0:.2f} %".format(conf * 100.0) @@ -4127,8 +3472,7 @@ class pol_map(object): self.an_int.remove() self.str_int = ( r"$F_{{\lambda}}^{{int}}$({0:.0f} $\AA$) = {1} $ergs \cdot cm^{{-2}} \cdot s^{{-1}} \cdot \AA^{{-1}}$".format( - self.pivot_wav, - sci_not(I_reg * self.map_convert, I_reg_err * self.map_convert, 2), + self.pivot_wav, sci_not(I_reg * self.map_convert, I_reg_err * self.map_convert, 2) ) + "\n" + r"$P^{{int}}$ = {0:.1f} $\pm$ {1:.1f} %".format(P_reg * 100.0, np.ceil(P_reg_err * 1000.0) / 10.0) @@ -4158,19 +3502,13 @@ class pol_map(object): horizontalalignment="left", ) if self.region is not None: - self.cont = ax.contour( - self.region.astype(float), - levels=[0.5], - colors="white", - linewidths=0.8, - ) + self.cont = ax.contour(self.region.astype(float), levels=[0.5], colors="white", linewidths=0.8) fig.canvas.draw_idle() return self.an_int else: str_int = ( r"$F_{{\lambda}}^{{int}}$({0:.0f} $\AA$) = {1} $ergs \cdot cm^{{-2}} \cdot s^{{-1}} \cdot \AA^{{-1}}$".format( - self.pivot_wav, - sci_not(I_reg * self.map_convert, I_reg_err * self.map_convert, 2), + self.pivot_wav, sci_not(I_reg * self.map_convert, I_reg_err * self.map_convert, 2) ) + "\n" + r"$P^{{int}}$ = {0:.1f} $\pm$ {1:.1f} %".format(P_reg * 100.0, np.ceil(P_reg_err * 1000.0) / 10.0) @@ -4200,12 +3538,7 @@ class pol_map(object): horizontalalignment="left", ) if self.region is not None: - ax.contour( - self.region.astype(float), - levels=[0.5], - colors="white", - linewidths=0.8, - ) + ax.contour(self.region.astype(float), levels=[0.5], colors="white", linewidths=0.8) fig.canvas.draw_idle() @@ -4213,84 +3546,21 @@ if __name__ == "__main__": import argparse parser = argparse.ArgumentParser(description="Interactively plot the pipeline products") + parser.add_argument("-f", "--file", metavar="path", required=False, help="The full or relative path to the data product", type=str, default=None) + parser.add_argument("-p", "--pcut", metavar="p_cut", required=False, help="The cut in signal-to-noise for the polarization degree", type=float, default=3.0) + parser.add_argument("-i", "--snri", metavar="snri_cut", required=False, help="The cut in signal-to-noise for the intensity", type=float, default=3.0) parser.add_argument( - "-f", - "--file", - metavar="path", - required=False, - help="The full or relative path to the data product", - type=str, - default=None, + "-st", "--step-vec", metavar="step_vec", required=False, help="Quantity of vectors to be shown, 1 is all, 2 is every other, etc.", type=int, default=1 ) parser.add_argument( - "-p", - "--pcut", - metavar="p_cut", - required=False, - help="The cut in signal-to-noise for the polarization degree", - type=float, - default=3.0, + "-sc", "--scale-vec", metavar="scale_vec", required=False, help="Size of the 100%% polarization vector in pixel units", type=float, default=3.0 ) + parser.add_argument("-pa", "--pang-err", action="store_true", required=False, help="Whether the polarization angle uncertainties should be displayed") + parser.add_argument("-l", "--lim", metavar="flux_lim", nargs=2, required=False, help="Limits for the intensity map", type=float, default=None) parser.add_argument( - "-i", - "--snri", - metavar="snri_cut", - required=False, - help="The cut in signal-to-noise for the intensity", - type=float, - default=3.0, - ) - parser.add_argument( - "-st", - "--step-vec", - metavar="step_vec", - required=False, - help="Quantity of vectors to be shown, 1 is all, 2 is every other, etc.", - type=int, - default=1, - ) - parser.add_argument( - "-sc", - "--scale-vec", - metavar="scale_vec", - required=False, - help="Size of the 100%% polarization vector in pixel units", - type=float, - default=3.0, - ) - parser.add_argument( - "-pa", - "--pang-err", - action="store_true", - required=False, - help="Whether the polarization angle uncertainties should be displayed", - ) - parser.add_argument( - "-l", - "--lim", - metavar="flux_lim", - nargs=2, - required=False, - help="Limits for the intensity map", - type=float, - default=None, - ) - parser.add_argument( - "-pdf", - "--static-pdf", - metavar="static_pdf", - required=False, - help="Whether the analysis tool or the static pdfs should be outputed", - default=None, - ) - parser.add_argument( - "-t", - "--type", - metavar="type", - required=False, - help="Type of plot to be be outputed", - default=None, + "-pdf", "--static-pdf", metavar="static_pdf", required=False, help="Whether the analysis tool or the static pdfs should be outputed", default=None ) + parser.add_argument("-t", "--type", metavar="type", required=False, help="Type of plot to be be outputed", default=None) args = parser.parse_args() if args.file is not None: diff --git a/package/lib/reduction.py b/package/lib/reduction.py index 65639fe..2898a5d 100755 --- a/package/lib/reduction.py +++ b/package/lib/reduction.py @@ -190,8 +190,8 @@ def bin_ndarray(ndarray, new_shape, operation="sum"): Example ------- - >>> m = np.arange(0,100,1).reshape((10,10)) - >>> n = bin_ndarray(m, new_shape=(5,5), operation='sum') + >>> m = np.arange(0, 100, 1).reshape((10, 10)) + >>> n = bin_ndarray(m, new_shape=(5, 5), operation="sum") >>> print(n) [[ 22 30 38 46 54] @@ -278,9 +278,7 @@ def crop_array(data_array, headers, error_array=None, data_mask=None, step=5, nu if null_val is None: null_val = [1.00 * error.mean() for error in error_array] elif type(null_val) is float: - null_val = [ - null_val, - ] * error_array.shape[0] + null_val = [null_val] * error_array.shape[0] vertex = np.zeros((data_array.shape[0], 4), dtype=int) for i, image in enumerate(data_array): # Get vertex of the rectangular convex hull of each image @@ -349,10 +347,7 @@ def crop_array(data_array, headers, error_array=None, data_mask=None, step=5, nu headers, vmin=convert_flux * data_array[data_array > 0.0].mean() / 5.0, vmax=convert_flux * data_array[data_array > 0.0].max(), - rectangle=[ - rectangle, - ] - * len(headers), + rectangle=[rectangle] * len(headers), savename=savename + "_crop_region", plots_folder=plots_folder, ) @@ -632,12 +627,7 @@ def rebin_array(data_array, error_array, headers, pxsize=2, scale="px", operatio # Compute binning ratio if scale.lower() in ["px", "pixel"]: - Dxy_arr[i] = np.array( - [ - pxsize, - ] - * 2 - ) + Dxy_arr[i] = np.array([pxsize] * 2) scale = "px" elif scale.lower() in ["arcsec", "arcseconds"]: Dxy_arr[i] = np.array(pxsize / np.abs(w.wcs.cdelt) / 3600.0) @@ -947,12 +937,7 @@ def smooth_data(data_array, error_array, data_mask, headers, FWHM=1.5, scale="pi dist_rc = np.where(data_mask, np.sqrt((r - xx) ** 2 + (c - yy) ** 2), fmax) # Catch expected "OverflowWarning" as we overflow values that are not in the image with warnings.catch_warnings(record=True) as w: - g_rc = np.array( - [ - np.exp(-0.5 * (dist_rc / stdev) ** 2) / (2.0 * np.pi * stdev**2), - ] - * data_array.shape[0] - ) + g_rc = np.array([np.exp(-0.5 * (dist_rc / stdev) ** 2) / (2.0 * np.pi * stdev**2)] * data_array.shape[0]) # Apply weighted combination smoothed[r, c] = np.where(data_mask[r, c], np.sum(data_array * weight * g_rc) / np.sum(weight * g_rc), data_array.mean(axis=0)[r, c]) error[r, c] = np.where( @@ -1447,9 +1432,7 @@ def compute_Stokes(data_array, error_array, data_mask, headers, FWHM=None, scale all_Q_stokes = np.zeros((np.unique(rotate).size, data_array.shape[1], data_array.shape[2])) all_U_stokes = np.zeros((np.unique(rotate).size, data_array.shape[1], data_array.shape[2])) all_Stokes_cov = np.zeros((np.unique(rotate).size, 3, 3, data_array.shape[1], data_array.shape[2])) - all_header_stokes = [ - {}, - ] * np.unique(rotate).size + all_header_stokes = [{}] * np.unique(rotate).size for i, rot in enumerate(np.unique(rotate)): rot_mask = rotate == rot From 70035a96264cde0e88c326235234b9c404daf599 Mon Sep 17 00:00:00 2001 From: Thibault Barnouin Date: Fri, 14 Mar 2025 14:30:30 +0100 Subject: [PATCH 18/33] Save the raw total flux image as PrimaryHDU fix rebase display on main --- package/FOC_reduction.py | 24 +++++++++++++++++------ package/lib/fits.py | 41 ++++++++++++++++++++++++++++++++++------ 2 files changed, 53 insertions(+), 12 deletions(-) diff --git a/package/FOC_reduction.py b/package/FOC_reduction.py index 55bd806..053706f 100755 --- a/package/FOC_reduction.py +++ b/package/FOC_reduction.py @@ -40,8 +40,8 @@ def main(target=None, proposal_id=None, infiles=None, output_dir="./data", crop= display_crop = False # Background estimation - error_sub_type = "scott" # sqrt, sturges, rice, scott, freedman-diaconis (default) or shape (example (51, 51)) - subtract_error = 2.0 + error_sub_type = "scott" # sqrt, sturges, rice, freedman-diaconis, scott (default) or shape (example (51, 51)) + subtract_error = 3.0 display_bkg = True # Data binning @@ -182,6 +182,14 @@ def main(target=None, proposal_id=None, infiles=None, output_dir="./data", crop= norm=LogNorm(vmin=data_array[data_array > 0.0].min() * headers[0]["photflam"], vmax=data_array[data_array > 0.0].max() * headers[0]["photflam"]), ) + flux_data, flux_error, flux_mask, flux_head = ( + deepcopy(data_array.sum(axis=0)), + deepcopy(np.sqrt(np.sum(error_array**2, axis=0))), + deepcopy(data_mask), + deepcopy(headers[0]), + ) + flux_head["EXPTIME"] = np.sum([head["EXPTIME"] for head in headers]) + # Rebin data to desired pixel size. if (pxsize is not None) and not (pxsize == 1 and pxscale.lower() in ["px", "pixel", "pixels"]): data_array, error_array, headers, Dxy, data_mask = proj_red.rebin_array( @@ -233,6 +241,8 @@ def main(target=None, proposal_id=None, infiles=None, output_dir="./data", crop= I_bkg, Q_bkg, U_bkg, S_cov_bkg, data_mask_bkg, header_bkg = proj_red.rotate_Stokes( I_bkg, Q_bkg, U_bkg, S_cov_bkg, np.array(True).reshape(1, 1), header_bkg, SNRi_cut=None ) + flux_data, flux_error, flux_mask, flux_head = proj_red.rotate_data(np.array([flux_data]), np.array([flux_error]), flux_mask, [flux_head]) + flux_data, flux_error, flux_head = flux_data[0], flux_error[0], flux_head[0] # Compute polarimetric parameters (polarization degree and angle). P, debiased_P, s_P, s_P_P, PA, s_PA, s_PA_P = proj_red.compute_pol(I_stokes, Q_stokes, U_stokes, Stokes_cov, header_stokes) @@ -258,8 +268,10 @@ def main(target=None, proposal_id=None, infiles=None, output_dir="./data", crop= figname, data_folder=data_folder, return_hdul=True, + flux_data=flux_data, + flux_head=flux_head, ) - outfiles.append("/".join([data_folder, Stokes_hdul[0].header["FILENAME"] + ".fits"])) + outfiles.append("/".join([data_folder, Stokes_hdul["I_STOKES"].header["FILENAME"] + ".fits"])) # Step 5: # crop to desired region of interest (roi) @@ -269,15 +281,15 @@ def main(target=None, proposal_id=None, infiles=None, output_dir="./data", crop= stokescrop.crop() stokescrop.write_to("/".join([data_folder, figname + ".fits"])) Stokes_hdul, header_stokes = stokescrop.hdul_crop, stokescrop.hdul_crop[0].header - outfiles.append("/".join([data_folder, Stokes_hdul[0].header["FILENAME"] + ".fits"])) + outfiles.append("/".join([data_folder, Stokes_hdul["I_STOKES"].header["FILENAME"] + ".fits"])) data_mask = Stokes_hdul["data_mask"].data.astype(bool) print( "F_int({0:.0f} Angs) = ({1} ± {2})e{3} ergs.cm^-2.s^-1.Angs^-1".format( header_stokes["PHOTPLAM"], *sci_not( - Stokes_hdul[0].data[data_mask].sum() * header_stokes["PHOTFLAM"], - np.sqrt(Stokes_hdul[3].data[0, 0][data_mask].sum()) * header_stokes["PHOTFLAM"], + Stokes_hdul["I_STOKES"].data[data_mask].sum() * header_stokes["PHOTFLAM"], + np.sqrt(Stokes_hdul["IQU_COV_MATRIX"].data[0, 0][data_mask].sum()) * header_stokes["PHOTFLAM"], 2, out=int, ), diff --git a/package/lib/fits.py b/package/lib/fits.py index 4e5aef7..99c41ed 100755 --- a/package/lib/fits.py +++ b/package/lib/fits.py @@ -101,7 +101,24 @@ def get_obs_data(infiles, data_folder="", compute_flux=False): def save_Stokes( - I_stokes, Q_stokes, U_stokes, Stokes_cov, P, debiased_P, s_P, s_P_P, PA, s_PA, s_PA_P, header_stokes, data_mask, filename, data_folder="", return_hdul=False + I_stokes, + Q_stokes, + U_stokes, + Stokes_cov, + P, + debiased_P, + s_P, + s_P_P, + PA, + s_PA, + s_PA_P, + header_stokes, + data_mask, + filename, + data_folder="", + return_hdul=False, + flux_data=None, + flux_head=None, ): """ Save computed polarimetry parameters to a single fits file, @@ -194,11 +211,23 @@ def save_Stokes( hdul = fits.HDUList([]) # Add I_stokes as PrimaryHDU - header["datatype"] = ("I_stokes", "type of data stored in the HDU") - I_stokes[(1 - data_mask).astype(bool)] = 0.0 - primary_hdu = fits.PrimaryHDU(data=I_stokes, header=header) - primary_hdu.name = "I_stokes" - hdul.append(primary_hdu) + if flux_data is None: + header["datatype"] = ("I_stokes", "type of data stored in the HDU") + I_stokes[(1 - data_mask).astype(bool)] = 0.0 + primary_hdu = fits.PrimaryHDU(data=I_stokes, header=header) + primary_hdu.name = "I_stokes" + hdul.append(primary_hdu) + else: + flux_head["TELESCOP"], flux_head["INSTRUME"] = header["TELESCOP"], header["INSTRUME"] + header["datatype"] = ("Flux map", "type of data stored in the HDU") + primary_hdu = fits.PrimaryHDU(data=flux_data, header=flux_head) + primary_hdu.name = "Flux map" + hdul.append(primary_hdu) + header["datatype"] = ("I_stokes", "type of data stored in the HDU") + I_stokes[(1 - data_mask).astype(bool)] = 0.0 + image_hdu = fits.ImageHDU(data=I_stokes, header=header) + image_hdu.name = "I_stokes" + hdul.append(image_hdu) # Add Q, U, Stokes_cov, P, s_P, PA, s_PA to the HDUList for data, name in [ From e3a3abdb0d5aec3f4cd9c0e26a78ad3cb4d4e3fc Mon Sep 17 00:00:00 2001 From: Thibault Barnouin Date: Wed, 19 Mar 2025 16:38:03 +0100 Subject: [PATCH 19/33] change histogram binning to numpy function --- package/lib/background.py | 32 +++++++++----------------------- 1 file changed, 9 insertions(+), 23 deletions(-) diff --git a/package/lib/background.py b/package/lib/background.py index 99a09f2..96509d2 100755 --- a/package/lib/background.py +++ b/package/lib/background.py @@ -278,7 +278,7 @@ def bkg_fit(data, error, mask, headers, subtract_error=True, display=False, save return n_data_array, n_error_array, headers, background -def bkg_hist(data, error, mask, headers, sub_type=None, subtract_error=True, display=False, savename=None, plots_folder=""): +def bkg_hist(data, error, mask, headers, n_bins=None, subtract_error=True, display=False, savename=None, plots_folder=""): """ ---------- Inputs: @@ -333,29 +333,15 @@ def bkg_hist(data, error, mask, headers, sub_type=None, subtract_error=True, dis for i, image in enumerate(data): # Compute the Count-rate histogram for the image n_mask = np.logical_and(mask, image > 0.0) - if sub_type is not None: - if isinstance(sub_type, int): - n_bins = sub_type - elif sub_type.lower() in ["square-root", "squareroot", "sqrt"]: - n_bins = np.fix(np.sqrt(image[n_mask].size)).astype(int) # Square-root - elif sub_type.lower() in ["sturges"]: - n_bins = np.ceil(np.log2(image[n_mask].size)).astype(int) + 1 # Sturges - elif sub_type.lower() in ["rice"]: - n_bins = 2 * np.fix(np.power(image[n_mask].size, 1 / 3)).astype(int) # Rice - elif sub_type.lower() in ["freedman-diaconis", "freedmandiaconis", "freedman", "diaconis"]: - n_bins = np.fix( - (image[n_mask].max() - image[n_mask].min()) - / (2 * np.subtract(*np.percentile(image[n_mask], [75, 25])) / np.power(image[n_mask].size, 1 / 3)) - ).astype(int) # Freedman-Diaconis - else: # Fallback - n_bins = np.fix((image[n_mask].max() - image[n_mask].min()) / (3.5 * image[n_mask].std() / np.power(image[n_mask].size, 1 / 3))).astype( - int - ) # Scott - else: # Default statistic - n_bins = np.fix((image[n_mask].max() - image[n_mask].min()) / (3.5 * image[n_mask].std() / np.power(image[n_mask].size, 1 / 3))).astype( - int - ) # Scott + if not isinstance(n_bins, int) and n_bins not in ["auto", "fd", "doane", "scott", "stone", "rice", "sturges", "sqrt"]: + match n_bins.lower(): + case "square-root" | "squareroot": + n_bins = "sqrt" + case "freedman-diaconis" | "freedmandiaconis": + n_bins = "fd" + case _: + n_bins = "scott" hist, bin_edges = np.histogram(np.log(image[n_mask]), bins=n_bins) histograms.append(hist) binning.append(np.exp(bin_centers(bin_edges))) From fb63b389e165757bb7e2ccf753f0047fb19c0543 Mon Sep 17 00:00:00 2001 From: Thibault Barnouin Date: Wed, 19 Mar 2025 16:38:36 +0100 Subject: [PATCH 20/33] save raw flux in fits file and display fix rebase display on main --- package/FOC_reduction.py | 37 ++++--- package/lib/fits.py | 26 +++-- package/lib/plots.py | 234 ++++++++++++++++++++++++++------------- package/lib/reduction.py | 17 ++- 4 files changed, 204 insertions(+), 110 deletions(-) diff --git a/package/FOC_reduction.py b/package/FOC_reduction.py index 053706f..899c9e6 100755 --- a/package/FOC_reduction.py +++ b/package/FOC_reduction.py @@ -20,6 +20,7 @@ import lib.reduction as proj_red # Functions used in reduction pipeline import numpy as np from lib.utils import princ_angle, sci_not from matplotlib.colors import LogNorm +from astropy.wcs import WCS def main(target=None, proposal_id=None, infiles=None, output_dir="./data", crop=False, interactive=False): @@ -40,9 +41,9 @@ def main(target=None, proposal_id=None, infiles=None, output_dir="./data", crop= display_crop = False # Background estimation - error_sub_type = "scott" # sqrt, sturges, rice, freedman-diaconis, scott (default) or shape (example (51, 51)) + error_sub_type = "scott" # auto, sqrt, sturges, rice, freedman-diaconis, scott (default) or shape (example (51, 51)) subtract_error = 3.0 - display_bkg = True + display_bkg = False # Data binning pxsize = 4 @@ -63,6 +64,7 @@ def main(target=None, proposal_id=None, infiles=None, output_dir="./data", crop= smoothing_scale = "px" # pixel or arcsec # Rotation + rotate_data = False rotate_North = True # Polarization map output @@ -142,7 +144,7 @@ def main(target=None, proposal_id=None, infiles=None, output_dir="./data", crop= ) # Rotate data to have same orientation - rotate_data = np.unique([np.round(float(head["ORIENTAT"]), 3) for head in headers]).size != 1 + rotate_data = rotate_data or np.unique([np.round(float(head["ORIENTAT"]), 3) for head in headers]).size != 1 if rotate_data: ang = np.mean([head["ORIENTAT"] for head in headers]) for head in headers: @@ -158,7 +160,6 @@ def main(target=None, proposal_id=None, infiles=None, output_dir="./data", crop= vmin=data_array[data_array > 0.0].min() * headers[0]["photflam"], vmax=data_array[data_array > 0.0].max() * headers[0]["photflam"] ), ) - # Align and rescale images with oversampling. data_array, error_array, headers, data_mask, shifts, error_shifts = proj_red.align_data( data_array, @@ -243,6 +244,8 @@ def main(target=None, proposal_id=None, infiles=None, output_dir="./data", crop= ) flux_data, flux_error, flux_mask, flux_head = proj_red.rotate_data(np.array([flux_data]), np.array([flux_error]), flux_mask, [flux_head]) flux_data, flux_error, flux_head = flux_data[0], flux_error[0], flux_head[0] + elif not rotate_data: + figname += "_noROT" # Compute polarimetric parameters (polarization degree and angle). P, debiased_P, s_P, s_P_P, PA, s_PA, s_PA_P = proj_red.compute_pol(I_stokes, Q_stokes, U_stokes, Stokes_cov, header_stokes) @@ -250,7 +253,7 @@ def main(target=None, proposal_id=None, infiles=None, output_dir="./data", crop= # Step 4: # Save image to FITS. - figname = "_".join([figname, figtype]) if figtype != "" else figname + savename = "_".join([figname, figtype]) if figtype != "" else figname Stokes_hdul = proj_fits.save_Stokes( I_stokes, Q_stokes, @@ -265,18 +268,18 @@ def main(target=None, proposal_id=None, infiles=None, output_dir="./data", crop= s_PA_P, header_stokes, data_mask, - figname, + savename, data_folder=data_folder, return_hdul=True, - flux_data=flux_data, + flux_data=np.array([flux_data, flux_error, flux_mask]), flux_head=flux_head, ) - outfiles.append("/".join([data_folder, Stokes_hdul["I_STOKES"].header["FILENAME"] + ".fits"])) + outfiles.append("/".join([data_folder, Stokes_hdul[0].header["FILENAME"] + ".fits"])) # Step 5: # crop to desired region of interest (roi) if crop: - figname += "_crop" + savename += "_crop" stokescrop = proj_plots.crop_Stokes(deepcopy(Stokes_hdul), norm=LogNorm()) stokescrop.crop() stokescrop.write_to("/".join([data_folder, figname + ".fits"])) @@ -286,10 +289,10 @@ def main(target=None, proposal_id=None, infiles=None, output_dir="./data", crop= data_mask = Stokes_hdul["data_mask"].data.astype(bool) print( "F_int({0:.0f} Angs) = ({1} ± {2})e{3} ergs.cm^-2.s^-1.Angs^-1".format( - header_stokes["PHOTPLAM"], + flux_head["PHOTPLAM"], *sci_not( - Stokes_hdul["I_STOKES"].data[data_mask].sum() * header_stokes["PHOTFLAM"], - np.sqrt(Stokes_hdul["IQU_COV_MATRIX"].data[0, 0][data_mask].sum()) * header_stokes["PHOTFLAM"], + flux_data[flux_mask].sum() * flux_head["PHOTFLAM"], + np.sqrt(np.sum(flux_error[flux_mask] ** 2)) * flux_head["PHOTFLAM"], 2, out=int, ), @@ -315,12 +318,12 @@ def main(target=None, proposal_id=None, infiles=None, output_dir="./data", crop= flux_lim=flux_lim, step_vec=step_vec, scale_vec=scale_vec, - savename="_".join([figname]), + savename=figname, plots_folder=plots_folder, ) for figtype, figsuffix in zip( - ["Intensity", "Pol_flux", "Pol_deg", "Pol_ang", "I_err", "P_err", "SNRi", "SNRp", "confp"], - ["I", "P_flux", "P", "PA", "I_err", "P_err", "SNRi", "SNRp", "confP"], + ["FluxDensity", "Intensity", "Pol_flux", "Pol_deg", "Pol_ang", "I_err", "P_err", "SNRi", "SNRp", "confp"], + ["F", "I", "P_flux", "P", "PA", "I_err", "P_err", "SNRi", "SNRp", "confP"], ): try: proj_plots.polarization_map( @@ -331,7 +334,7 @@ def main(target=None, proposal_id=None, infiles=None, output_dir="./data", crop= flux_lim=flux_lim, step_vec=step_vec, scale_vec=scale_vec, - savename="_".join([figname, figsuffix]), + savename="_".join([savename, figsuffix]), plots_folder=plots_folder, display=figtype, ) @@ -339,7 +342,7 @@ def main(target=None, proposal_id=None, infiles=None, output_dir="./data", crop= pass elif not interactive: proj_plots.polarization_map( - deepcopy(Stokes_hdul), data_mask, P_cut=P_cut, SNRi_cut=SNRi_cut, savename=figname, plots_folder=plots_folder, display="integrate" + deepcopy(Stokes_hdul), data_mask, P_cut=P_cut, SNRi_cut=SNRi_cut, savename=savename, plots_folder=plots_folder, display="integrate" ) elif pxscale.lower() not in ["full", "integrate"]: proj_plots.pol_map(Stokes_hdul, P_cut=P_cut, SNRi_cut=SNRi_cut, step_vec=step_vec, scale_vec=scale_vec, flux_lim=flux_lim) diff --git a/package/lib/fits.py b/package/lib/fits.py index 99c41ed..67a62d9 100755 --- a/package/lib/fits.py +++ b/package/lib/fits.py @@ -47,9 +47,9 @@ def get_obs_data(infiles, data_folder="", compute_flux=False): wcs_array.append(WCS(header=f[0].header, fobj=f).celestial) f.flush() # Save pixel area for flux density computation - if headers[i]["PXFORMT"] == "NORMAL": + if "PXFORMT" in headers[i].keys() and headers[i]["PXFORMT"] == "NORMAL": headers[i]["PXAREA"] = 1.96e-4 # 14x14 milliarcsec squared pixel area in arcsec^2 - elif headers[i]["PXFORMT"] == "ZOOM": + elif "PXFORMT" in headers[i].keys() and headers[i]["PXFORMT"] == "ZOOM": headers[i]["PXAREA"] = 4.06e-4 # 29x14 milliarcsec squared pixel area in arcsec^2 else: headers[i]["PXAREA"] = 1.0 # unknown default to 1 arcsec^2 @@ -85,10 +85,10 @@ def get_obs_data(infiles, data_folder="", compute_flux=False): # force WCS for POL60 to have same pixel size as POL0 and POL120 is_pol60 = np.array([head["filtnam1"].lower() == "pol60" for head in headers], dtype=bool) cdelt = np.round(np.array([WCS(head).wcs.cdelt[:2] for head in headers]), 10) - if np.unique(cdelt[np.logical_not(is_pol60)], axis=0).size != 2: + if np.any(is_pol60) and np.unique(cdelt[np.logical_not(is_pol60)], axis=0).size != 2: print(np.unique(cdelt[np.logical_not(is_pol60)], axis=0)) raise ValueError("Not all images have same pixel size") - else: + elif np.any(is_pol60): for i in np.arange(len(headers))[is_pol60]: headers[i]["cdelt1"], headers[i]["cdelt2"] = np.unique(cdelt[np.logical_not(is_pol60)], axis=0)[0] @@ -174,7 +174,7 @@ def save_Stokes( header["PROPOSID"] = (header_stokes["PROPOSID"], "PEP proposal identifier for observation") header["TARGNAME"] = (header_stokes["TARGNAME"], "Target name") header["ORIENTAT"] = (header_stokes["ORIENTAT"], "Angle between North and the y-axis of the image") - header["FILENAME"] = (filename, "ORIGINAL FILENAME") + header["FILENAME"] = (filename, "Original filename") header["BKG_TYPE"] = (header_stokes["BKG_TYPE"], "Bkg estimation method used during reduction") header["BKG_SUB"] = (header_stokes["BKG_SUB"], "Amount of bkg subtracted from images") header["SMOOTH"] = (header_stokes["SMOOTH"] if "SMOOTH" in list(header_stokes.keys()) else "None", "Smoothing method used during reduction") @@ -210,7 +210,7 @@ def save_Stokes( # Create HDUList object hdul = fits.HDUList([]) - # Add I_stokes as PrimaryHDU + # Add Flux density as PrimaryHDU if flux_data is None: header["datatype"] = ("I_stokes", "type of data stored in the HDU") I_stokes[(1 - data_mask).astype(bool)] = 0.0 @@ -218,10 +218,16 @@ def save_Stokes( primary_hdu.name = "I_stokes" hdul.append(primary_hdu) else: - flux_head["TELESCOP"], flux_head["INSTRUME"] = header["TELESCOP"], header["INSTRUME"] - header["datatype"] = ("Flux map", "type of data stored in the HDU") - primary_hdu = fits.PrimaryHDU(data=flux_data, header=flux_head) - primary_hdu.name = "Flux map" + flux_head["FILENAME"] = header["FILENAME"] + head = WCS(flux_head).deepcopy().to_header() + for key in [key for key in header.keys() if key not in ["SMOOTH", "SAMPLING"]]: + try: + head[key] = flux_head[key] + except KeyError: + head[key] = header[key] + header["datatype"] = ("Flux_density", "type of data stored in the HDU") + primary_hdu = fits.PrimaryHDU(data=flux_data, header=head) + primary_hdu.name = "Flux_density" hdul.append(primary_hdu) header["datatype"] = ("I_stokes", "type of data stored in the HDU") I_stokes[(1 - data_mask).astype(bool)] = 0.0 diff --git a/package/lib/plots.py b/package/lib/plots.py index c5a92d2..f5dbbd6 100755 --- a/package/lib/plots.py +++ b/package/lib/plots.py @@ -133,6 +133,12 @@ def plot_obs(data_array, headers, rectangle=None, shifts=None, savename=None, pl ax_curr.arrow(x, y, dx, dy, length_includes_head=True, width=0.1, head_width=0.3, color="g") ax_curr.plot([x, x], [0, data.shape[0] - 1], "--", lw=2, color="g", alpha=0.85) ax_curr.plot([0, data.shape[1] - 1], [y, y], "--", lw=2, color="g", alpha=0.85) + # position of centroid + ax_curr.plot([data.shape[1] / 2, data.shape[1] / 2], [0, data.shape[0] - 1], "--", lw=2, color="b", alpha=0.85) + ax_curr.plot([0, data.shape[1] - 1], [data.shape[0] / 2, data.shape[0] / 2], "--", lw=2, color="b", alpha=0.85) + cr_x, cr_y = head["CRPIX1"], head["CRPIX2"] + # Plot WCS reference point + ax_curr.plot([cr_x], [cr_y], "+", lw=2, color="r", alpha=0.85) if rectangle is not None: x, y, width, height, angle, color = rectangle[i] ax_curr.add_patch(Rectangle((x, y), width, height, angle=angle, edgecolor=color, fill=False)) @@ -189,7 +195,7 @@ def plot_Stokes(Stokes, savename=None, plots_folder=""): for dataset in [stkI, stkQ, stkU]: dataset[np.logical_not(data_mask)] = np.nan - wcs = WCS(Stokes[0]).deepcopy() + wcs = WCS(Stokes["I_STOKES"]).deepcopy() # Plot figure plt.rcParams.update({"font.size": 14}) @@ -288,6 +294,9 @@ def polarization_map( The figure and ax created for interactive contour maps. """ # Get data + flux = Stokes[0].data[0].copy() * Stokes[0].header["PHOTFLAM"] + flux_error = Stokes[0].data[1].copy() * Stokes[0].header["PHOTFLAM"] + flux_mask = Stokes[0].data[2].astype(bool).copy() stkI = Stokes["I_stokes"].data.copy() stkQ = Stokes["Q_stokes"].data.copy() stkU = Stokes["U_stokes"].data.copy() @@ -302,6 +311,20 @@ def polarization_map( data_mask = np.zeros(stkI.shape).astype(bool) data_mask[stkI > 0.0] = True + wcs = WCS(Stokes["I_STOKES"]).deepcopy() + pivot_wav = Stokes["I_STOKES"].header["photplam"] + convert_flux = Stokes["I_STOKES"].header["photflam"] + + # Get integrated flux values from sum + I_diluted = stkI[data_mask].sum() * convert_flux + I_diluted_err = np.sqrt(np.sum(stk_cov[0, 0][data_mask])) * convert_flux + + # Get integrated polarization values from header + P_diluted = Stokes["I_STOKES"].header["P_int"] + P_diluted_err = Stokes["I_STOKES"].header["sP_int"] + PA_diluted = Stokes["I_STOKES"].header["PA_int"] + PA_diluted_err = Stokes["I_STOKES"].header["sPA_int"] + # Compute confidence level map QN, UN, QN_ERR, UN_ERR = np.full((4, stkI.shape[0], stkI.shape[1]), np.nan) for nflux, sflux in zip([QN, UN, QN_ERR, UN_ERR], [stkQ, stkU, np.sqrt(stk_cov[1, 1]), np.sqrt(stk_cov[2, 2])]): @@ -314,12 +337,8 @@ def polarization_map( for j in range(3): stk_cov[i][j][np.logical_not(data_mask)] = np.nan - wcs = WCS(Stokes[0]).deepcopy() - pivot_wav = Stokes[0].header["photplam"] - convert_flux = Stokes[0].header["photflam"] - # Plot Stokes parameters map - if display is None or display.lower() in ["default"]: + if display is None or display.lower() in ["pol", "polarization", "polarisation", "pol_deg", "p"]: plot_Stokes(Stokes, savename=savename, plots_folder=plots_folder) # Compute SNR and apply cuts @@ -363,7 +382,7 @@ def polarization_map( fig = plt.figure(figsize=(7 * ratiox, 7 * ratioy), layout="constrained") if ax is None: ax = fig.add_subplot(111, projection=wcs) - ax.set(aspect="equal", fc="k") + ax.set(aspect="equal", fc="k", xlim=(0, stkI.shape[1]), ylim=(0, stkI.shape[0])) # fig.subplots_adjust(hspace=0, wspace=0, left=0.102, right=1.02) # ax.coords.grid(True, color='white', ls='dotted', alpha=0.5) @@ -409,7 +428,25 @@ def polarization_map( ax.set_facecolor("white") font_color = "black" - if display.lower() in ["i", "intensity"]: + if display.lower() in ["f", "flux", "fluxdensity"]: + # If no display selected, show intensity map + display = "f" + if flux_lim is not None: + vmin, vmax = flux_lim + else: + vmin, vmax = np.max(flux[flux > 0.0]) / 2e3, np.max(flux[flux > 0.0]) + imflux, cr = flux.copy(), WCS(Stokes[0].header).wcs.crpix.astype(int) + imflux[cr[1] - 1 : cr[1] + 2, cr[0] - 1 : cr[0] + 2] = np.nan + im = ax.imshow( + imflux, transform=ax.get_transform(WCS(Stokes[0].header).celestial), norm=LogNorm(vmin, vmax), aspect="equal", cmap=kwargs["cmap"], alpha=1.0 + ) + fig.colorbar(im, ax=ax, aspect=50, shrink=0.60, pad=0.025, label=r"$F_{\lambda}$ [$ergs \cdot cm^{-2} \cdot s^{-1} \cdot \AA$]") + levelsF = np.array([0.8, 2.0, 5.0, 10.0, 20.0, 50.0]) / 100.0 * vmax + print("Flux density contour levels : ", levelsF) + ax.contour(flux, levels=levelsF, transform=ax.get_transform(WCS(Stokes[0].header).celestial), colors="grey", linewidths=0.5) + ax.plot(*WCS(Stokes[1]).wcs.crpix, "g+") + I_diluted, I_diluted_err = np.sum(flux[flux_mask]), np.sqrt(np.sum(flux_error[flux_mask] ** 2)) + elif display.lower() in ["i", "intensity"]: # If no display selected, show intensity map display = "i" if flux_lim is not None: @@ -420,9 +457,12 @@ def polarization_map( vmin, vmax = 1.0 / 2.0 * np.median(np.sqrt(stk_cov[0, 0][stkI > 0.0]) * convert_flux), np.max(stkI[stkI > 0.0] * convert_flux) im = ax.imshow(stkI * convert_flux, norm=LogNorm(vmin, vmax), aspect="equal", cmap=kwargs["cmap"], alpha=1.0) fig.colorbar(im, ax=ax, aspect=50, shrink=0.60, pad=0.025, label=r"$F_{\lambda}$ [$ergs \cdot cm^{-2} \cdot s^{-1} \cdot \AA^{-1}$]") - levelsI = np.array([0.8, 2.0, 5.0, 10.0, 20.0, 50.0]) / 100.0 * vmax - print("Flux density contour levels : ", levelsI) - ax.contour(stkI * convert_flux, levels=levelsI, colors="grey", linewidths=0.5) + # levelsI = np.array([0.8, 2.0, 5.0, 10.0, 20.0, 50.0]) / 100.0 * vmax + # print("Stokes I contour levels : ", levelsI) + # ax.contour(stkI * convert_flux, levels=levelsI, colors="grey", linewidths=0.5) + levelsF = np.array([0.8, 2.0, 5.0, 10.0, 20.0, 50.0]) / 100.0 * np.max(flux[flux > 0.0]) + print("Flux density contour levels : ", levelsF) + ax.contour(flux, levels=levelsF, transform=ax.get_transform(WCS(Stokes[0].header).celestial), colors="grey", linewidths=0.5) elif display.lower() in ["pf", "pol_flux"]: # Display polarization flux display = "pf" @@ -512,22 +552,18 @@ def polarization_map( # Defaults to intensity map if flux_lim is not None: vmin, vmax = flux_lim - elif mask.sum() > 0.0: - vmin, vmax = 1.0 * np.mean(np.sqrt(stk_cov[0, 0][mask]) * convert_flux), np.max(stkI[stkI > 0.0] * convert_flux) else: - vmin, vmax = 1.0 * np.mean(np.sqrt(stk_cov[0, 0][stkI > 0.0]) * convert_flux), np.max(stkI[stkI > 0.0] * convert_flux) - im = ax.imshow(stkI * convert_flux, norm=LogNorm(vmin, vmax), aspect="equal", cmap=kwargs["cmap"], alpha=1.0) + vmin, vmax = np.max(flux[flux > 0.0] * convert_flux) / 2e3, np.max(flux[flux > 0.0] * convert_flux) + im = ax.imshow( + flux * Stokes[0].header["PHOTFLAM"], + transform=ax.get_transform(WCS(Stokes[0].header).celestial), + norm=LogNorm(vmin, vmax), + aspect="equal", + cmap=kwargs["cmap"], + alpha=1.0, + ) fig.colorbar(im, ax=ax, aspect=50, shrink=0.60, pad=0.025, label=r"$F_{\lambda}$ [$ergs \cdot cm^{-2} \cdot s^{-1} \cdot \AA$]") - - # Get integrated flux values from sum - I_diluted = stkI[data_mask].sum() - I_diluted_err = np.sqrt(np.sum(stk_cov[0, 0][data_mask])) - - # Get integrated polarization values from header - P_diluted = Stokes[0].header["P_int"] - P_diluted_err = Stokes[0].header["sP_int"] - PA_diluted = Stokes[0].header["PA_int"] - PA_diluted_err = Stokes[0].header["sPA_int"] + I_diluted, I_diluted_err = np.sum(flux[flux_mask]), np.sqrt(np.sum(flux_error[flux_mask] ** 2)) plt.rcParams.update({"font.size": 11}) px_size = wcs.wcs.get_cdelt()[0] * 3600.0 @@ -545,12 +581,12 @@ def polarization_map( back_length=0.0, head_length=7.0, head_width=7.0, - angle=-Stokes[0].header["orientat"], + angle=-Stokes["I_STOKES"].header["orientat"], text_props={"ec": "k", "fc": font_color, "alpha": 1, "lw": 0.5}, arrow_props={"ec": "k", "fc": font_color, "alpha": 1, "lw": 1}, ) - if display.lower() in ["i", "s_i", "snri", "pf", "p", "pa", "s_p", "snrp", "confp"] and step_vec != 0: + if display.lower() in ["f", "i", "s_i", "snri", "pf", "p", "pa", "s_p", "snrp", "confp"] and step_vec != 0: if scale_vec == -1: poldata[np.isfinite(poldata)] = 1.0 / 2.0 step_vec = 1 @@ -1174,6 +1210,8 @@ class overplot_chandra(align_maps): other_data = deepcopy(self.other_data) other_wcs = self.other_wcs.deepcopy() if zoom != 1: + from scipy.ndimage import zoom as sc_zoom + other_data = sc_zoom(other_data, zoom) other_wcs.wcs.crpix *= zoom other_wcs.wcs.cdelt /= zoom @@ -2439,9 +2477,9 @@ class pol_map(object): self.conf = PCconf(self.QN, self.UN, self.QN_ERR, self.UN_ERR) # Get data - self.targ = self.Stokes[0].header["targname"] - self.pivot_wav = self.Stokes[0].header["photplam"] - self.map_convert = self.Stokes[0].header["photflam"] + self.targ = self.Stokes["I_STOKES"].header["targname"] + self.pivot_wav = self.Stokes["I_STOKES"].header["photplam"] + self.map_convert = self.Stokes["I_STOKES"].header["photflam"] # Create figure plt.rcParams.update({"font.size": 10}) @@ -2535,7 +2573,7 @@ class pol_map(object): def select_roi(event): if self.data is None: - self.data = self.Stokes[0].data + self.data = self.Stokes["I_STOKES"].data if self.selected: self.selected = False self.region = deepcopy(self.select_instance.mask.astype(bool)) @@ -2577,7 +2615,7 @@ class pol_map(object): def select_aperture(event): if self.data is None: - self.data = self.Stokes[0].data + self.data = self.Stokes["I_STOKES"].data if self.selected: self.selected = False self.select_instance.update_mask() @@ -2634,7 +2672,7 @@ class pol_map(object): def select_slit(event): if self.data is None: - self.data = self.Stokes[0].data + self.data = self.Stokes["I_STOKES"].data if self.selected: self.selected = False self.select_instance.update_mask() @@ -2911,7 +2949,19 @@ class pol_map(object): @property def wcs(self): - return WCS(self.Stokes[0].header).celestial.deepcopy() + return WCS(self.Stokes["I_STOKES"].header).celestial.deepcopy() + + @property + def Flux(self): + return self.Stokes[0].data[0] * self.Stokes[0].header["PHOTFLAM"] + + @property + def Flux_err(self): + return self.Stokes[0].data[1] * self.Stokes[0].header["PHOTFLAM"] + + @property + def Flux_mask(self): + return self.Stokes[0].data[2].astype(bool) @property def I(self): @@ -2975,7 +3025,7 @@ class pol_map(object): @property def data_mask(self): - return self.Stokes["DATA_MASK"].data + return self.Stokes["DATA_MASK"].data.astype(bool) def set_data_mask(self, mask): self.Stokes[np.argmax([self.Stokes[i].header["datatype"] == "Data_mask" for i in range(len(self.Stokes))])].data = mask.astype(float) @@ -3046,7 +3096,7 @@ class pol_map(object): back_length=0.0, head_length=10.0, head_width=10.0, - angle=-self.Stokes[0].header["orientat"], + angle=-self.Stokes["I_STOKES"].header["orientat"], color="white", text_props={"ec": None, "fc": "w", "alpha": 1, "lw": 0.4}, arrow_props={"ec": None, "fc": "w", "alpha": 1, "lw": 1}, @@ -3054,18 +3104,20 @@ class pol_map(object): ax.add_artist(self.north_dir) def display(self, fig=None, ax=None, flux_lim=None): - norm = None - if self.display_selection is None: - self.display_selection = "total_flux" + kwargs = dict([]) if flux_lim is None: flux_lim = self.flux_lim - if self.display_selection.lower() in ["total_flux"]: - self.data = self.I * self.map_convert + if self.display_selection is None or self.display_selection.lower() in ["total_flux"]: + self.data = self.Flux # self.I * self.map_convert if flux_lim is None: vmin, vmax = (1.0 / 2.0 * np.median(self.data[self.data > 0.0]), np.max(self.data[self.data > 0.0])) else: vmin, vmax = flux_lim - norm = LogNorm(vmin, vmax) + kwargs["norm"] = LogNorm(vmin, vmax) + if ax is None: + kwargs["transform"] = self.ax.get_transform(WCS(self.Stokes[0].header).celestial) + else: + kwargs["transform"] = ax.get_transform(WCS(self.Stokes[0].header).celestial) label = r"$F_{\lambda}$ [$ergs \cdot cm^{-2} \cdot s^{-1} \cdot \AA^{-1}$]" elif self.display_selection.lower() in ["pol_flux"]: self.data = self.I * self.map_convert * self.P @@ -3073,28 +3125,28 @@ class pol_map(object): vmin, vmax = (1.0 / 2.0 * np.median(self.I[self.I > 0.0] * self.map_convert), np.max(self.I[self.I > 0.0] * self.map_convert)) else: vmin, vmax = flux_lim - norm = LogNorm(vmin, vmax) + kwargs["norm"] = LogNorm(vmin, vmax) label = r"$P \cdot F_{\lambda}$ [$ergs \cdot cm^{-2} \cdot s^{-1} \cdot \AA^{-1}$]" elif self.display_selection.lower() in ["pol_deg"]: self.data = self.P * 100.0 - vmin, vmax = 0.0, np.max(self.data[self.P > self.s_P]) + kwargs["vmin"], kwargs["vmax"] = 0.0, np.max(self.data[self.P > self.s_P]) label = r"$P$ [%]" elif self.display_selection.lower() in ["pol_ang"]: self.data = princ_angle(self.PA) - vmin, vmax = 0, 180.0 + kwargs["vmin"], kwargs["vmax"] = 0, 180.0 label = r"$\theta_{P}$ [°]" elif self.display_selection.lower() in ["snri"]: s_I = np.sqrt(self.IQU_cov[0, 0]) SNRi = np.zeros(self.I.shape) SNRi[s_I > 0.0] = self.I[s_I > 0.0] / s_I[s_I > 0.0] self.data = SNRi - vmin, vmax = 0.0, np.max(self.data[self.data > 0.0]) + kwargs["vmin"], kwargs["vmax"] = 0.0, np.max(self.data[self.data > 0.0]) label = r"$I_{Stokes}/\sigma_{I}$" elif self.display_selection.lower() in ["snrp"]: SNRp = np.zeros(self.P.shape) SNRp[self.s_P > 0.0] = self.P[self.s_P > 0.0] / self.s_P[self.s_P > 0.0] self.data = SNRp - vmin, vmax = 0.0, np.max(self.data[self.data > 0.0]) + kwargs["vmin"], kwargs["vmax"] = 0.0, np.max(self.data[self.data > 0.0]) label = r"$P/\sigma_{P}$" if fig is None: @@ -3105,22 +3157,17 @@ class pol_map(object): self.cbar.remove() if hasattr(self, "im"): self.im.remove() - if norm is not None: - self.im = ax.imshow(self.data, norm=norm, aspect="equal", cmap="inferno") - else: - self.im = ax.imshow(self.data, vmin=vmin, vmax=vmax, aspect="equal", cmap="inferno") + self.im = ax.imshow(self.data, aspect="equal", cmap="inferno", **kwargs) + ax.set(xlim=(0, self.I.shape[1]), ylim=(0, self.I.shape[0])) plt.rcParams.update({"font.size": 14}) self.cbar = fig.colorbar(self.im, ax=ax, aspect=50, shrink=0.75, pad=0.025, label=label) plt.rcParams.update({"font.size": 10}) fig.canvas.draw_idle() return self.im else: - if norm is not None: - im = ax.imshow(self.data, norm=norm, aspect="equal", cmap="inferno") - else: - im = ax.imshow(self.data, vmin=vmin, vmax=vmax, aspect="equal", cmap="inferno") - ax.set_xlim(0, self.data.shape[1]) - ax.set_ylim(0, self.data.shape[0]) + im = ax.imshow(self.data, aspect="equal", cmap="inferno", **kwargs) + ax.set_xlim(0, self.I.shape[1]) + ax.set_ylim(0, self.I.shape[0]) plt.rcParams.update({"font.size": 14}) fig.colorbar(im, ax=ax, aspect=50, shrink=0.75, pad=0.025, label=label) plt.rcParams.update({"font.size": 10}) @@ -3144,24 +3191,24 @@ class pol_map(object): ax = self.ax if hasattr(self, "quiver"): self.quiver.remove() - self.quiver = ax.quiver( - X[:: self.step_vec, :: self.step_vec], - Y[:: self.step_vec, :: self.step_vec], - XY_U[:: self.step_vec, :: self.step_vec], - XY_V[:: self.step_vec, :: self.step_vec], - units="xy", - scale=1.0 / scale_vec, - scale_units="xy", - pivot="mid", - headwidth=0.0, - headlength=0.0, - headaxislength=0.0, - width=0.3, - linewidth=0.6, - color="white", - edgecolor="black", - ) if self.pa_err: + self.quiver = ax.quiver( + X[:: self.step_vec, :: self.step_vec], + Y[:: self.step_vec, :: self.step_vec], + XY_U[:: self.step_vec, :: self.step_vec], + XY_V[:: self.step_vec, :: self.step_vec], + units="xy", + scale=1.0 / scale_vec, + scale_units="xy", + pivot="mid", + headwidth=0.0, + headlength=0.0, + headaxislength=0.0, + width=0.1, + # linewidth=0.6, + color="black", + edgecolor="black", + ) XY_U_err1, XY_V_err1 = ( P_cut * np.cos(np.pi / 2.0 + (self.PA + 3.0 * self.s_PA) * np.pi / 180.0), P_cut * np.sin(np.pi / 2.0 + (self.PA + 3.0 * self.s_PA) * np.pi / 180.0), @@ -3210,6 +3257,25 @@ class pol_map(object): edgecolor="black", ls="dashed", ) + else: + self.quiver = ax.quiver( + X[:: self.step_vec, :: self.step_vec], + Y[:: self.step_vec, :: self.step_vec], + XY_U[:: self.step_vec, :: self.step_vec], + XY_V[:: self.step_vec, :: self.step_vec], + units="xy", + scale=1.0 / scale_vec, + scale_units="xy", + pivot="mid", + headwidth=0.0, + headlength=0.0, + headaxislength=0.0, + width=0.3, + linewidth=0.6, + color="white", + edgecolor="black", + ) + fig.canvas.draw_idle() return self.quiver else: @@ -3281,12 +3347,20 @@ class pol_map(object): str_conf = "" if self.region is None: s_I = np.sqrt(self.IQU_cov[0, 0]) - I_reg = self.I.sum() - I_reg_err = np.sqrt(np.sum(s_I**2)) - P_reg = self.Stokes[0].header["P_int"] - P_reg_err = self.Stokes[0].header["sP_int"] - PA_reg = self.Stokes[0].header["PA_int"] - PA_reg_err = self.Stokes[0].header["sPA_int"] + I_reg = ( + np.sum(self.Flux[self.Flux_mask]) / self.map_convert + if self.display_selection is None or self.display_selection.lower() in ["total_flux"] + else np.sum(self.I[self.data_mask]) + ) + I_reg_err = ( + np.sqrt(np.sum(self.Flux_err[self.Flux_mask] ** 2)) / self.map_convert + if self.display_selection is None or self.display_selection.lower() in ["total_flux"] + else np.sqrt(np.sum(s_I[self.data_mask] ** 2)) + ) + P_reg = self.Stokes["I_STOKES"].header["P_int"] + P_reg_err = self.Stokes["I_STOKES"].header["sP_int"] + PA_reg = self.Stokes["I_STOKES"].header["PA_int"] + PA_reg_err = self.Stokes["I_STOKES"].header["sPA_int"] s_I = np.sqrt(self.IQU_cov[0, 0]) s_Q = np.sqrt(self.IQU_cov[1, 1]) diff --git a/package/lib/reduction.py b/package/lib/reduction.py index 89d52b5..0f4e5d2 100755 --- a/package/lib/reduction.py +++ b/package/lib/reduction.py @@ -52,7 +52,6 @@ from scipy.ndimage import rotate as sc_rotate from scipy.ndimage import shift as sc_shift from scipy.signal import fftconvolve -from .background import bkg_fit, bkg_hist, bkg_mini from .convex_hull import image_hull from .cross_correlation import phase_cross_correlation from .deconvolve import deconvolve_im, gaussian2d, gaussian_psf, zeropad @@ -308,6 +307,8 @@ def crop_array(data_array, headers, error_array=None, data_mask=None, step=5, nu # Update CRPIX value in the associated header curr_wcs = WCS(crop_headers[i]).celestial.deepcopy() curr_wcs.wcs.crpix[:2] = curr_wcs.wcs.crpix[:2] - np.array([v_array[2], v_array[0]]) + curr_wcs.array_shape = crop_array[i].shape + curr_wcs.wcs.set() crop_headers[i].update(curr_wcs.to_header()) crop_headers[i]["naxis1"], crop_headers[i]["naxis2"] = crop_array[i].shape @@ -523,20 +524,22 @@ def get_error( # estimated to less than 3% err_flat = data * 0.03 + from .background import bkg_fit, bkg_hist, bkg_mini + if sub_type is None: n_data_array, c_error_bkg, headers, background = bkg_hist( data, error, mask, headers, subtract_error=subtract_error, display=display, savename=savename, plots_folder=plots_folder ) sub_type, subtract_error = "histogram ", str(int(subtract_error > 0.0)) elif isinstance(sub_type, str): - if sub_type.lower() in ["auto"]: + if sub_type.lower() in ["fit"]: n_data_array, c_error_bkg, headers, background = bkg_fit( data, error, mask, headers, subtract_error=subtract_error, display=display, savename=savename, plots_folder=plots_folder ) sub_type, subtract_error = "histogram fit ", "mean+%.1fsigma" % subtract_error if subtract_error != 0.0 else 0.0 else: n_data_array, c_error_bkg, headers, background = bkg_hist( - data, error, mask, headers, sub_type=sub_type, subtract_error=subtract_error, display=display, savename=savename, plots_folder=plots_folder + data, error, mask, headers, n_bins=sub_type, subtract_error=subtract_error, display=display, savename=savename, plots_folder=plots_folder ) sub_type, subtract_error = "histogram ", "mean+%.1fsigma" % subtract_error if subtract_error != 0.0 else 0.0 elif isinstance(sub_type, tuple): @@ -665,6 +668,7 @@ def rebin_array(data_array, error_array, headers, pxsize=2, scale="px", operatio nw.wcs.cdelt *= Dxy nw.wcs.crpix /= Dxy nw.array_shape = new_shape + nw.wcs.set() new_header["NAXIS1"], new_header["NAXIS2"] = nw.array_shape new_header["PXAREA"] *= Dxy[0] * Dxy[1] for key, val in nw.to_header().items(): @@ -844,7 +848,10 @@ def align_data( new_crpix = np.array([wcs.wcs.crpix for wcs in headers_wcs]) + shifts[:, ::-1] + res_shift[::-1] for i in range(len(headers_wcs)): headers_wcs[i].wcs.crpix = new_crpix[0] + headers_wcs[i].array_shape = (res_shape, res_shape) + headers_wcs[i].wcs.set() headers[i].update(headers_wcs[i].to_header()) + headers[i]["NAXIS1"], headers[i]["NAXIS2"] = res_shape, res_shape data_mask = rescaled_mask.all(axis=0) data_array, error_array, data_mask, headers = crop_array(rescaled_image, headers, rescaled_error, data_mask, null_val=0.01 * background) @@ -1706,9 +1713,11 @@ def rotate_Stokes(I_stokes, Q_stokes, U_stokes, Stokes_cov, data_mask, header_st new_wcs.wcs.pc = np.dot(mrot, new_wcs.wcs.pc) new_wcs.wcs.crpix = np.dot(mrot, new_wcs.wcs.crpix - old_center[::-1]) + new_center[::-1] + new_wcs.array_shape = shape new_wcs.wcs.set() for key, val in new_wcs.to_header().items(): new_header_stokes.set(key, val) + new_header_stokes["NAXIS1"], new_header_stokes["NAXIS2"] = new_wcs.array_shape new_header_stokes["ORIENTAT"] += ang # Nan handling : @@ -1809,9 +1818,11 @@ def rotate_data(data_array, error_array, data_mask, headers): new_wcs = WCS(header).celestial.deepcopy() new_wcs.wcs.pc[:2, :2] = np.dot(mrot, new_wcs.wcs.pc[:2, :2]) new_wcs.wcs.crpix[:2] = np.dot(mrot, new_wcs.wcs.crpix[:2] - old_center[::-1]) + new_center[::-1] + new_wcs.array_shape = shape new_wcs.wcs.set() for key, val in new_wcs.to_header().items(): new_header[key] = val + new_header["NAXIS1"], new_header["NAXIS2"] = new_wcs.array_shape new_header["ORIENTAT"] = np.arccos(new_wcs.celestial.wcs.pc[0, 0]) * 180.0 / np.pi new_header["ROTATE"] = ang new_headers.append(new_header) From 6457b4ffe1a8437d016325293f70eaa9a8edfce0 Mon Sep 17 00:00:00 2001 From: Thibault Barnouin Date: Tue, 1 Apr 2025 16:20:40 +0200 Subject: [PATCH 21/33] looking for displacement of WCS in pipeline --- package/FOC_reduction.py | 1 - package/lib/plots.py | 2 +- package/lib/reduction.py | 4 +++- 3 files changed, 4 insertions(+), 3 deletions(-) diff --git a/package/FOC_reduction.py b/package/FOC_reduction.py index 899c9e6..b42f57f 100755 --- a/package/FOC_reduction.py +++ b/package/FOC_reduction.py @@ -20,7 +20,6 @@ import lib.reduction as proj_red # Functions used in reduction pipeline import numpy as np from lib.utils import princ_angle, sci_not from matplotlib.colors import LogNorm -from astropy.wcs import WCS def main(target=None, proposal_id=None, infiles=None, output_dir="./data", crop=False, interactive=False): diff --git a/package/lib/plots.py b/package/lib/plots.py index f5dbbd6..21bf79c 100755 --- a/package/lib/plots.py +++ b/package/lib/plots.py @@ -456,12 +456,12 @@ def polarization_map( else: vmin, vmax = 1.0 / 2.0 * np.median(np.sqrt(stk_cov[0, 0][stkI > 0.0]) * convert_flux), np.max(stkI[stkI > 0.0] * convert_flux) im = ax.imshow(stkI * convert_flux, norm=LogNorm(vmin, vmax), aspect="equal", cmap=kwargs["cmap"], alpha=1.0) + ax.plot(*WCS(Stokes[1]).wcs.crpix, "g+") fig.colorbar(im, ax=ax, aspect=50, shrink=0.60, pad=0.025, label=r"$F_{\lambda}$ [$ergs \cdot cm^{-2} \cdot s^{-1} \cdot \AA^{-1}$]") # levelsI = np.array([0.8, 2.0, 5.0, 10.0, 20.0, 50.0]) / 100.0 * vmax # print("Stokes I contour levels : ", levelsI) # ax.contour(stkI * convert_flux, levels=levelsI, colors="grey", linewidths=0.5) levelsF = np.array([0.8, 2.0, 5.0, 10.0, 20.0, 50.0]) / 100.0 * np.max(flux[flux > 0.0]) - print("Flux density contour levels : ", levelsF) ax.contour(flux, levels=levelsF, transform=ax.get_transform(WCS(Stokes[0].header).celestial), colors="grey", linewidths=0.5) elif display.lower() in ["pf", "pol_flux"]: # Display polarization flux diff --git a/package/lib/reduction.py b/package/lib/reduction.py index 0f4e5d2..2898a5d 100755 --- a/package/lib/reduction.py +++ b/package/lib/reduction.py @@ -637,7 +637,8 @@ def rebin_array(data_array, error_array, headers, pxsize=2, scale="px", operatio pxsize, scale = "", "full" else: raise ValueError("'{0:s}' invalid scale for binning.".format(scale)) - new_shape = np.ceil(min(image.shape / Dxy_arr, key=lambda x: x[0] + x[1])).astype(int) + new_shape_float = min(image.shape / Dxy_arr, key=lambda x: x[0] + x[1]) + new_shape = np.ceil(new_shape_float).astype(int) for i, (image, error, header) in enumerate(list(zip(data_array, error_array, headers))): # Get current pixel size @@ -666,6 +667,7 @@ def rebin_array(data_array, error_array, headers, pxsize=2, scale="px", operatio # Update header nw = w.deepcopy() nw.wcs.cdelt *= Dxy + # nw.wcs.crpix += np.abs(new_shape_float - new_shape) * np.array(new_shape) / Dxy nw.wcs.crpix /= Dxy nw.array_shape = new_shape nw.wcs.set() From e1409167b985f01a78e7a74082704839d074a2cc Mon Sep 17 00:00:00 2001 From: Thibault Barnouin Date: Tue, 11 Mar 2025 16:14:14 +0100 Subject: [PATCH 22/33] small improvments to ConfCenter et emission_center --- package/lib/utils.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/package/lib/utils.py b/package/lib/utils.py index ee9fc5b..b76b063 100755 --- a/package/lib/utils.py +++ b/package/lib/utils.py @@ -109,7 +109,7 @@ def CenterConf(mask, PA, sPA): conf[np.isfinite(PA)] = 0.5 * np.exp(-0.5 * chi2[np.isfinite(PA)]) c0 = np.unravel_index(np.argmax(conf), conf.shape)[::-1] - result = minimize(chisq, c0, bounds=[(0, PA.shape[1]), (0.0, PA.shape[0])], method="trust-constr") + result = minimize(chisq, c0, bounds=[(0, PA.shape[1]), (0.0, PA.shape[0])], method="TNC") if result.success: print("Center of emission found: reduced chi_squared {0:.2f}/{1:d}={2:.2f}".format(chisq(result.x) * (Nobs - 2), Nobs - 2, chisq(result.x))) else: From 8bfaee25fb4a96014e560403c5738bacef41f721 Mon Sep 17 00:00:00 2001 From: Thibault Barnouin Date: Fri, 14 Mar 2025 14:20:38 +0100 Subject: [PATCH 23/33] rollback CenterConf to simple CDF of chi2 for 2 dof --- package/lib/utils.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/package/lib/utils.py b/package/lib/utils.py index b76b063..ee9fc5b 100755 --- a/package/lib/utils.py +++ b/package/lib/utils.py @@ -109,7 +109,7 @@ def CenterConf(mask, PA, sPA): conf[np.isfinite(PA)] = 0.5 * np.exp(-0.5 * chi2[np.isfinite(PA)]) c0 = np.unravel_index(np.argmax(conf), conf.shape)[::-1] - result = minimize(chisq, c0, bounds=[(0, PA.shape[1]), (0.0, PA.shape[0])], method="TNC") + result = minimize(chisq, c0, bounds=[(0, PA.shape[1]), (0.0, PA.shape[0])], method="trust-constr") if result.success: print("Center of emission found: reduced chi_squared {0:.2f}/{1:d}={2:.2f}".format(chisq(result.x) * (Nobs - 2), Nobs - 2, chisq(result.x))) else: From 1193b2b091997c2ea4c089bda7743ec88d70137c Mon Sep 17 00:00:00 2001 From: Thibault Barnouin Date: Tue, 1 Apr 2025 17:02:31 +0200 Subject: [PATCH 24/33] fix WCS computation, cdelt should not be sorted fix rebase display on main --- package/lib/fits.py | 2 +- package/lib/utils.py | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/package/lib/fits.py b/package/lib/fits.py index 67a62d9..cbacc85 100755 --- a/package/lib/fits.py +++ b/package/lib/fits.py @@ -16,7 +16,7 @@ from astropy.io import fits from astropy.wcs import WCS from .convex_hull import clean_ROI -from .utils import wcs_PA +from .utils import wcs_PA, princ_angle def get_obs_data(infiles, data_folder="", compute_flux=False): diff --git a/package/lib/utils.py b/package/lib/utils.py index ee9fc5b..9e36dbc 100755 --- a/package/lib/utils.py +++ b/package/lib/utils.py @@ -177,4 +177,4 @@ def wcs_PA(PC21, PC22): orient = np.arccos(PC22) * 180.0 / np.pi elif (abs(PC21) < abs(PC22)) and (PC22 < 0): orient = -np.arccos(PC22) * 180.0 / np.pi - return orient + return princ_angle(orient) From e597d81e4ff708ecc0607ce6b1dad1d649506eaa Mon Sep 17 00:00:00 2001 From: Thibault Barnouin Date: Fri, 14 Mar 2025 14:30:30 +0100 Subject: [PATCH 25/33] Save the raw total flux image as PrimaryHDU fast forward fixes to WCS fix forward merging of display --- package/FOC_reduction.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/package/FOC_reduction.py b/package/FOC_reduction.py index b42f57f..dcea9fa 100755 --- a/package/FOC_reduction.py +++ b/package/FOC_reduction.py @@ -273,7 +273,7 @@ def main(target=None, proposal_id=None, infiles=None, output_dir="./data", crop= flux_data=np.array([flux_data, flux_error, flux_mask]), flux_head=flux_head, ) - outfiles.append("/".join([data_folder, Stokes_hdul[0].header["FILENAME"] + ".fits"])) + outfiles.append("/".join([data_folder, Stokes_hdul["I_STOKES"].header["FILENAME"] + ".fits"])) # Step 5: # crop to desired region of interest (roi) From 4b104d79d44a5b8c9dd39556eafc4bdede54939e Mon Sep 17 00:00:00 2001 From: Thibault Barnouin Date: Wed, 19 Mar 2025 16:38:36 +0100 Subject: [PATCH 26/33] save raw flux in fits file and display forward fixes to display with WCS --- package/FOC_reduction.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/package/FOC_reduction.py b/package/FOC_reduction.py index dcea9fa..899c9e6 100755 --- a/package/FOC_reduction.py +++ b/package/FOC_reduction.py @@ -20,6 +20,7 @@ import lib.reduction as proj_red # Functions used in reduction pipeline import numpy as np from lib.utils import princ_angle, sci_not from matplotlib.colors import LogNorm +from astropy.wcs import WCS def main(target=None, proposal_id=None, infiles=None, output_dir="./data", crop=False, interactive=False): @@ -273,7 +274,7 @@ def main(target=None, proposal_id=None, infiles=None, output_dir="./data", crop= flux_data=np.array([flux_data, flux_error, flux_mask]), flux_head=flux_head, ) - outfiles.append("/".join([data_folder, Stokes_hdul["I_STOKES"].header["FILENAME"] + ".fits"])) + outfiles.append("/".join([data_folder, Stokes_hdul[0].header["FILENAME"] + ".fits"])) # Step 5: # crop to desired region of interest (roi) From cebf4df8a4c6cb189b8c5dd196fc06e42db5ad30 Mon Sep 17 00:00:00 2001 From: Thibault Barnouin Date: Mon, 7 Apr 2025 17:20:00 +0200 Subject: [PATCH 27/33] modifications to add raw image display --- package/FOC_reduction.py | 10 +++++----- package/__init__.py | 2 +- package/lib/background.py | 3 ++- package/lib/plots.py | 24 +++++++++++++----------- 4 files changed, 21 insertions(+), 18 deletions(-) diff --git a/package/FOC_reduction.py b/package/FOC_reduction.py index 899c9e6..314e2ff 100755 --- a/package/FOC_reduction.py +++ b/package/FOC_reduction.py @@ -46,8 +46,8 @@ def main(target=None, proposal_id=None, infiles=None, output_dir="./data", crop= display_bkg = False # Data binning - pxsize = 4 - pxscale = "px" # pixel, arcsec or full + pxsize = 0.10 + pxscale = "arcsec" # pixel, arcsec or full rebin_operation = "sum" # sum or average # Alignement @@ -60,8 +60,8 @@ def main(target=None, proposal_id=None, infiles=None, output_dir="./data", crop= # Smoothing smoothing_function = "combine" # gaussian_after, weighted_gaussian_after, gaussian, weighted_gaussian or combine - smoothing_FWHM = 1.5 # If None, no smoothing is done - smoothing_scale = "px" # pixel or arcsec + smoothing_FWHM = 0.150 # If None, no smoothing is done + smoothing_scale = "arcsec" # pixel or arcsec # Rotation rotate_data = False @@ -71,7 +71,7 @@ def main(target=None, proposal_id=None, infiles=None, output_dir="./data", crop= P_cut = 3 # if >=1.0 cut on the signal-to-noise else cut on the confidence level in Q, U SNRi_cut = 1.0 # I measurments with SNR>30, which implies an uncertainty in P of 4.7%. flux_lim = None # lowest and highest flux displayed on plot, defaults to bkg and maximum in cut if None - scale_vec = 3 + scale_vec = 5 step_vec = 1 # plot all vectors in the array. if step_vec = 2, then every other vector will be plotted if step_vec = 0 then all vectors are displayed at full length # Pipeline start diff --git a/package/__init__.py b/package/__init__.py index a8ca1bd..59ef14c 100644 --- a/package/__init__.py +++ b/package/__init__.py @@ -1,2 +1,2 @@ from .lib import * -from . import FOC_reduction +from .FOC_reduction import main diff --git a/package/lib/background.py b/package/lib/background.py index 96509d2..3e8405f 100755 --- a/package/lib/background.py +++ b/package/lib/background.py @@ -18,7 +18,6 @@ import matplotlib.dates as mdates import matplotlib.pyplot as plt import numpy as np from astropy.time import Time -from lib.plots import plot_obs from matplotlib.colors import LogNorm from matplotlib.patches import Rectangle from scipy.optimize import curve_fit @@ -136,6 +135,8 @@ def display_bkg(data, background, std_bkg, headers, histograms=None, binning=Non fig2.subplots_adjust(hspace=0, wspace=0, right=1.0) fig2.colorbar(im2, ax=ax2, location="right", shrink=0.60, aspect=50, pad=0.025, label=r"Flux [$ergs \cdot cm^{-2} \cdot s^{-1} \cdot \AA^{-1}$]") + from .plots import plot_obs + if savename is not None: this_savename = deepcopy(savename) if savename[-4:] not in [".png", ".jpg", ".pdf"]: diff --git a/package/lib/plots.py b/package/lib/plots.py index 21bf79c..3089564 100755 --- a/package/lib/plots.py +++ b/package/lib/plots.py @@ -440,7 +440,7 @@ def polarization_map( im = ax.imshow( imflux, transform=ax.get_transform(WCS(Stokes[0].header).celestial), norm=LogNorm(vmin, vmax), aspect="equal", cmap=kwargs["cmap"], alpha=1.0 ) - fig.colorbar(im, ax=ax, aspect=50, shrink=0.60, pad=0.025, label=r"$F_{\lambda}$ [$ergs \cdot cm^{-2} \cdot s^{-1} \cdot \AA$]") + fig.colorbar(im, ax=ax, aspect=50, shrink=0.60, pad=0.015, fraction=0.03, label=r"$F_{\lambda}$ [$ergs \cdot cm^{-2} \cdot s^{-1} \cdot \AA$]") levelsF = np.array([0.8, 2.0, 5.0, 10.0, 20.0, 50.0]) / 100.0 * vmax print("Flux density contour levels : ", levelsF) ax.contour(flux, levels=levelsF, transform=ax.get_transform(WCS(Stokes[0].header).celestial), colors="grey", linewidths=0.5) @@ -457,7 +457,7 @@ def polarization_map( vmin, vmax = 1.0 / 2.0 * np.median(np.sqrt(stk_cov[0, 0][stkI > 0.0]) * convert_flux), np.max(stkI[stkI > 0.0] * convert_flux) im = ax.imshow(stkI * convert_flux, norm=LogNorm(vmin, vmax), aspect="equal", cmap=kwargs["cmap"], alpha=1.0) ax.plot(*WCS(Stokes[1]).wcs.crpix, "g+") - fig.colorbar(im, ax=ax, aspect=50, shrink=0.60, pad=0.025, label=r"$F_{\lambda}$ [$ergs \cdot cm^{-2} \cdot s^{-1} \cdot \AA^{-1}$]") + fig.colorbar(im, ax=ax, aspect=50, shrink=0.60, pad=0.015, fraction=0.03, label=r"$F_{\lambda}$ [$ergs \cdot cm^{-2} \cdot s^{-1} \cdot \AA^{-1}$]") # levelsI = np.array([0.8, 2.0, 5.0, 10.0, 20.0, 50.0]) / 100.0 * vmax # print("Stokes I contour levels : ", levelsI) # ax.contour(stkI * convert_flux, levels=levelsI, colors="grey", linewidths=0.5) @@ -476,7 +476,9 @@ def polarization_map( vmin, vmax = 1.0 / 2.0 * np.median(np.sqrt(stk_cov[0, 0][stkI > 0.0]) * convert_flux), np.max(stkI[stkI > 0.0] * convert_flux) pfmax = (stkI[stkI > 0.0] * pol[stkI > 0.0] * convert_flux).max() im = ax.imshow(stkI * convert_flux * pol, norm=LogNorm(vmin, vmax), aspect="equal", cmap=kwargs["cmap"], alpha=1.0) - fig.colorbar(im, ax=ax, aspect=50, shrink=0.60, pad=0.025, label=r"$F_{\lambda} \cdot P$ [$ergs \cdot cm^{-2} \cdot s^{-1} \cdot \AA^{-1}$]") + fig.colorbar( + im, ax=ax, aspect=50, shrink=0.60, pad=0.015, fraction=0.03, label=r"$F_{\lambda} \cdot P$ [$ergs \cdot cm^{-2} \cdot s^{-1} \cdot \AA^{-1}$]" + ) # levelsPf = np.linspace(0.0.60, 0.50, 5) * pfmax levelsPf = np.array([1.73, 13.0, 33.0, 66.0]) / 100.0 * pfmax print("Polarized flux density contour levels : ", levelsPf) @@ -486,13 +488,13 @@ def polarization_map( display = "p" vmin, vmax = 0.0, 100.0 im = ax.imshow(pol * 100.0, vmin=vmin, vmax=vmax, aspect="equal", cmap=kwargs["cmap"], alpha=1.0) - fig.colorbar(im, ax=ax, aspect=50, shrink=0.60, pad=0.025, label=r"$P$ [%]") + fig.colorbar(im, ax=ax, aspect=50, shrink=0.60, pad=0.015, fraction=0.03, label=r"$P$ [%]") elif display.lower() in ["pa", "pang", "pol_ang"]: # Display polarization degree map display = "pa" vmin, vmax = 0.0, 180.0 im = ax.imshow(princ_angle(pang), vmin=vmin, vmax=vmax, aspect="equal", cmap=kwargs["cmap"], alpha=1.0) - fig.colorbar(im, ax=ax, aspect=50, shrink=0.60, pad=0.025, label=r"$\theta_P$ [°]") + fig.colorbar(im, ax=ax, aspect=50, shrink=0.60, pad=0.015, fraction=0.03, label=r"$\theta_P$ [°]") elif display.lower() in ["s_p", "pol_err", "pol_deg_err"]: # Display polarization degree error map display = "s_p" @@ -502,7 +504,7 @@ def polarization_map( else: vmin, vmax = 0.0, 100.0 im = ax.imshow(pol_err * 100.0, vmin=vmin, vmax=vmax, aspect="equal", cmap="inferno_r", alpha=1.0) - fig.colorbar(im, ax=ax, aspect=50, shrink=0.60, pad=0.025, label=r"$\sigma_P$ [%]") + fig.colorbar(im, ax=ax, aspect=50, shrink=0.60, pad=0.015, fraction=0.03, label=r"$\sigma_P$ [%]") elif display.lower() in ["s_i", "i_err"]: # Display intensity error map display = "s_i" @@ -514,7 +516,7 @@ def polarization_map( im = ax.imshow(np.sqrt(stk_cov[0, 0]) * convert_flux, norm=LogNorm(vmin, vmax), aspect="equal", cmap="inferno_r", alpha=1.0) else: im = ax.imshow(np.sqrt(stk_cov[0, 0]) * convert_flux, aspect="equal", cmap=kwargs["cmap"], alpha=1.0) - fig.colorbar(im, ax=ax, aspect=50, shrink=0.60, pad=0.025, label=r"$\sigma_I$ [$ergs \cdot cm^{-2} \cdot s^{-1} \cdot \AA^{-1}$]") + fig.colorbar(im, ax=ax, aspect=50, shrink=0.60, pad=0.015, fraction=0.03, label=r"$\sigma_I$ [$ergs \cdot cm^{-2} \cdot s^{-1} \cdot \AA^{-1}$]") elif display.lower() in ["snri"]: # Display I_stokes signal-to-noise map display = "snri" @@ -526,7 +528,7 @@ def polarization_map( ax.contour(SNRi, levels=levelsSNRi, colors="grey", linewidths=0.5) else: im = ax.imshow(SNRi, aspect="equal", cmap=kwargs["cmap"], alpha=1.0) - fig.colorbar(im, ax=ax, aspect=50, shrink=0.60, pad=0.025, label=r"$I_{Stokes}/\sigma_{I}$") + fig.colorbar(im, ax=ax, aspect=50, shrink=0.60, pad=0.015, fraction=0.03, label=r"$I_{Stokes}/\sigma_{I}$") elif display.lower() in ["snr", "snrp"]: # Display polarization degree signal-to-noise map display = "snrp" @@ -538,7 +540,7 @@ def polarization_map( ax.contour(SNRp, levels=levelsSNRp, colors="grey", linewidths=0.5) else: im = ax.imshow(SNRp, aspect="equal", cmap=kwargs["cmap"], alpha=1.0) - fig.colorbar(im, ax=ax, aspect=50, shrink=0.60, pad=0.025, label=r"$P/\sigma_{P}$") + fig.colorbar(im, ax=ax, aspect=50, shrink=0.60, pad=0.015, fraction=0.03, label=r"$P/\sigma_{P}$") elif display.lower() in ["conf", "confp"]: # Display polarization degree signal-to-noise map display = "confp" @@ -547,7 +549,7 @@ def polarization_map( levelsconfp = np.array([0.500, 0.900, 0.990, 0.999]) print("confp contour levels : ", levelsconfp) ax.contour(confP, levels=levelsconfp, colors="grey", linewidths=0.5) - fig.colorbar(im, ax=ax, aspect=50, shrink=0.60, pad=0.025, label=r"$Conf_{P}$") + fig.colorbar(im, ax=ax, aspect=50, shrink=0.60, pad=0.015, fraction=0.03, label=r"$Conf_{P}$") else: # Defaults to intensity map if flux_lim is not None: @@ -562,7 +564,7 @@ def polarization_map( cmap=kwargs["cmap"], alpha=1.0, ) - fig.colorbar(im, ax=ax, aspect=50, shrink=0.60, pad=0.025, label=r"$F_{\lambda}$ [$ergs \cdot cm^{-2} \cdot s^{-1} \cdot \AA$]") + fig.colorbar(im, ax=ax, aspect=50, shrink=0.60, pad=0.015, fraction=0.03, label=r"$F_{\lambda}$ [$ergs \cdot cm^{-2} \cdot s^{-1} \cdot \AA$]") I_diluted, I_diluted_err = np.sum(flux[flux_mask]), np.sqrt(np.sum(flux_error[flux_mask] ** 2)) plt.rcParams.update({"font.size": 11}) From 2f4dda688feb04590961d76e578668f710040c91 Mon Sep 17 00:00:00 2001 From: Thibault Barnouin Date: Tue, 11 Mar 2025 16:14:14 +0100 Subject: [PATCH 28/33] small improvments to ConfCenter et emission_center --- package/lib/utils.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/package/lib/utils.py b/package/lib/utils.py index 9e36dbc..e56a29e 100755 --- a/package/lib/utils.py +++ b/package/lib/utils.py @@ -109,7 +109,7 @@ def CenterConf(mask, PA, sPA): conf[np.isfinite(PA)] = 0.5 * np.exp(-0.5 * chi2[np.isfinite(PA)]) c0 = np.unravel_index(np.argmax(conf), conf.shape)[::-1] - result = minimize(chisq, c0, bounds=[(0, PA.shape[1]), (0.0, PA.shape[0])], method="trust-constr") + result = minimize(chisq, c0, bounds=[(0, PA.shape[1]), (0.0, PA.shape[0])], method="TNC") if result.success: print("Center of emission found: reduced chi_squared {0:.2f}/{1:d}={2:.2f}".format(chisq(result.x) * (Nobs - 2), Nobs - 2, chisq(result.x))) else: From 8a7429103031eb3d0d5516d2bc587602289c9850 Mon Sep 17 00:00:00 2001 From: Thibault Barnouin Date: Fri, 14 Mar 2025 14:20:38 +0100 Subject: [PATCH 29/33] rollback CenterConf to simple CDF of chi2 for 2 dof --- package/lib/utils.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/package/lib/utils.py b/package/lib/utils.py index e56a29e..9e36dbc 100755 --- a/package/lib/utils.py +++ b/package/lib/utils.py @@ -109,7 +109,7 @@ def CenterConf(mask, PA, sPA): conf[np.isfinite(PA)] = 0.5 * np.exp(-0.5 * chi2[np.isfinite(PA)]) c0 = np.unravel_index(np.argmax(conf), conf.shape)[::-1] - result = minimize(chisq, c0, bounds=[(0, PA.shape[1]), (0.0, PA.shape[0])], method="TNC") + result = minimize(chisq, c0, bounds=[(0, PA.shape[1]), (0.0, PA.shape[0])], method="trust-constr") if result.success: print("Center of emission found: reduced chi_squared {0:.2f}/{1:d}={2:.2f}".format(chisq(result.x) * (Nobs - 2), Nobs - 2, chisq(result.x))) else: From d5d153e8554ef3417d18ab9f6d48401d04686019 Mon Sep 17 00:00:00 2001 From: Thibault Barnouin Date: Tue, 1 Apr 2025 17:02:31 +0200 Subject: [PATCH 30/33] fix WCS computation, cdelt should not be sorted fix rebase display on main --- package/FOC_reduction.py | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/package/FOC_reduction.py b/package/FOC_reduction.py index 314e2ff..d3aee2c 100755 --- a/package/FOC_reduction.py +++ b/package/FOC_reduction.py @@ -46,8 +46,8 @@ def main(target=None, proposal_id=None, infiles=None, output_dir="./data", crop= display_bkg = False # Data binning - pxsize = 0.10 - pxscale = "arcsec" # pixel, arcsec or full + pxsize = 4 + pxscale = "px" # pixel, arcsec or full rebin_operation = "sum" # sum or average # Alignement @@ -60,8 +60,8 @@ def main(target=None, proposal_id=None, infiles=None, output_dir="./data", crop= # Smoothing smoothing_function = "combine" # gaussian_after, weighted_gaussian_after, gaussian, weighted_gaussian or combine - smoothing_FWHM = 0.150 # If None, no smoothing is done - smoothing_scale = "arcsec" # pixel or arcsec + smoothing_FWHM = 1.5 # If None, no smoothing is done + smoothing_scale = "px" # pixel or arcsec # Rotation rotate_data = False From a7a6cbad1ccfc3f301763d4c48d18f26995f7bc6 Mon Sep 17 00:00:00 2001 From: Thibault Barnouin Date: Fri, 14 Mar 2025 14:30:30 +0100 Subject: [PATCH 31/33] Save the raw total flux image as PrimaryHDU fix for rebase display on main fix rebase display on main --- package/FOC_reduction.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/package/FOC_reduction.py b/package/FOC_reduction.py index d3aee2c..5abef77 100755 --- a/package/FOC_reduction.py +++ b/package/FOC_reduction.py @@ -274,7 +274,7 @@ def main(target=None, proposal_id=None, infiles=None, output_dir="./data", crop= flux_data=np.array([flux_data, flux_error, flux_mask]), flux_head=flux_head, ) - outfiles.append("/".join([data_folder, Stokes_hdul[0].header["FILENAME"] + ".fits"])) + outfiles.append("/".join([data_folder, Stokes_hdul["I_STOKES"].header["FILENAME"] + ".fits"])) # Step 5: # crop to desired region of interest (roi) From 77ce42bdd412292b5632aafb8631a56af3743d2b Mon Sep 17 00:00:00 2001 From: Thibault Barnouin Date: Wed, 19 Mar 2025 16:38:36 +0100 Subject: [PATCH 32/33] save raw flux in fits file and display fix rebase display on main fix rebase display on main --- package/FOC_reduction.py | 2 +- package/lib/plots.py | 4 ++-- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/package/FOC_reduction.py b/package/FOC_reduction.py index 5abef77..d3aee2c 100755 --- a/package/FOC_reduction.py +++ b/package/FOC_reduction.py @@ -274,7 +274,7 @@ def main(target=None, proposal_id=None, infiles=None, output_dir="./data", crop= flux_data=np.array([flux_data, flux_error, flux_mask]), flux_head=flux_head, ) - outfiles.append("/".join([data_folder, Stokes_hdul["I_STOKES"].header["FILENAME"] + ".fits"])) + outfiles.append("/".join([data_folder, Stokes_hdul[0].header["FILENAME"] + ".fits"])) # Step 5: # crop to desired region of interest (roi) diff --git a/package/lib/plots.py b/package/lib/plots.py index 3089564..8f63280 100755 --- a/package/lib/plots.py +++ b/package/lib/plots.py @@ -2951,7 +2951,7 @@ class pol_map(object): @property def wcs(self): - return WCS(self.Stokes["I_STOKES"].header).celestial.deepcopy() + return WCS(self.Stokes[0].header).celestial.deepcopy() @property def Flux(self): @@ -3098,7 +3098,7 @@ class pol_map(object): back_length=0.0, head_length=10.0, head_width=10.0, - angle=-self.Stokes["I_STOKES"].header["orientat"], + angle=-self.Stokes[0].header["orientat"], color="white", text_props={"ec": None, "fc": "w", "alpha": 1, "lw": 0.4}, arrow_props={"ec": None, "fc": "w", "alpha": 1, "lw": 1}, From fcbd5b8aec897e9711be00be3ab5d2efa174fe1d Mon Sep 17 00:00:00 2001 From: Thibault Barnouin Date: Tue, 8 Apr 2025 19:48:46 +0200 Subject: [PATCH 33/33] some more formatting fix rebase display on main --- package/FOC_reduction.py | 7 +------ 1 file changed, 1 insertion(+), 6 deletions(-) diff --git a/package/FOC_reduction.py b/package/FOC_reduction.py index d3aee2c..d3782cb 100755 --- a/package/FOC_reduction.py +++ b/package/FOC_reduction.py @@ -290,12 +290,7 @@ def main(target=None, proposal_id=None, infiles=None, output_dir="./data", crop= print( "F_int({0:.0f} Angs) = ({1} ± {2})e{3} ergs.cm^-2.s^-1.Angs^-1".format( flux_head["PHOTPLAM"], - *sci_not( - flux_data[flux_mask].sum() * flux_head["PHOTFLAM"], - np.sqrt(np.sum(flux_error[flux_mask] ** 2)) * flux_head["PHOTFLAM"], - 2, - out=int, - ), + *sci_not(flux_data[flux_mask].sum() * flux_head["PHOTFLAM"], np.sqrt(np.sum(flux_error[flux_mask] ** 2)) * flux_head["PHOTFLAM"], 2, out=int), ) ) print("P_int = {0:.1f} ± {1:.1f} %".format(header_stokes["p_int"] * 100.0, np.ceil(header_stokes["sP_int"] * 1000.0) / 10.0))