diff --git a/package/lib/plots.py b/package/lib/plots.py index dffee43..2ee4fbc 100755 --- a/package/lib/plots.py +++ b/package/lib/plots.py @@ -41,24 +41,29 @@ prototypes : """ from copy import deepcopy -import numpy as np from os.path import join as path_join -import matplotlib.pyplot as plt -from matplotlib.patches import Rectangle, Circle, FancyArrowPatch -from matplotlib.path import Path -from matplotlib.widgets import RectangleSelector, LassoSelector, Button, Slider, TextBox -from matplotlib.colors import LogNorm + import matplotlib.font_manager as fm import matplotlib.patheffects as pe -from mpl_toolkits.axes_grid1.anchored_artists import AnchoredSizeBar, AnchoredDirectionArrows -from astropy.wcs import WCS -from astropy.io import fits +import matplotlib.pyplot as plt +import numpy as np from astropy.coordinates import SkyCoord +from astropy.io import fits +from astropy.wcs import WCS +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 scipy.ndimage import zoom as sc_zoom + try: - from .utils import rot2D, princ_angle, sci_not + from .utils import princ_angle, rot2D, sci_not except ImportError: - from utils import rot2D, princ_angle, sci_not + from utils import princ_angle, rot2D, sci_not def plot_obs(data_array, headers, rectangle=None, savename=None, plots_folder="", **kwargs): @@ -91,19 +96,18 @@ def plot_obs(data_array, headers, rectangle=None, savename=None, plots_folder="" be saved. Not used if savename is None. Defaults to current folder. """ - 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']]) + 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]), dpi=200, layout='constrained', - sharex=True, sharey=True) + fig, ax = plt.subplots(shape[0], shape[1], figsize=(3 * shape[1], 3 * shape[0]), dpi=200, layout="constrained", sharex=True, sharey=True) r_pol = dict(pol0=0, pol60=1, pol120=2) c_pol = dict(pol0=0, pol60=0, pol120=0) for i, (data, head) in enumerate(zip(data_array, headers)): - instr = head['instrume'] - rootname = head['rootname'] - exptime = head['exptime'] - filt = head['filtnam1'] - convert = head['photflam'] + instr = head["instrume"] + rootname = head["rootname"] + exptime = head["exptime"] + filt = head["filtnam1"] + convert = head["photflam"] r_ax, c_ax = r_pol[filt.lower()], c_pol[filt.lower()] c_pol[filt.lower()] += 1 if shape[1] != 1: @@ -111,11 +115,11 @@ def plot_obs(data_array, headers, rectangle=None, savename=None, plots_folder="" else: ax_curr = ax[r_ax] # plots - if ('vmin' in kwargs.keys() or 'vmax' in kwargs.keys()): - vmin, vmax = kwargs['vmin'], kwargs['vmax'] - del kwargs['vmin'], kwargs['vmax'] + if "vmin" in kwargs.keys() or "vmax" in kwargs.keys(): + vmin, vmax = kwargs["vmin"], kwargs["vmax"] + del kwargs["vmin"], kwargs["vmax"] else: - vmin, vmax = convert*data[data > 0.].min()/10., convert*data[data > 0.].max() + vmin, vmax = convert * data[data > 0.0].min() / 10.0, convert * data[data > 0.0].max() for key, value in [["cmap", [["cmap", "gray"]]], ["norm", [["norm", LogNorm(vmin, vmax)]]]]: try: _ = kwargs[key] @@ -123,26 +127,28 @@ def plot_obs(data_array, headers, rectangle=None, savename=None, plots_folder="" for key_i, val_i in value: kwargs[key_i] = val_i # im = ax[r_ax][c_ax].imshow(convert*data, origin='lower', **kwargs) - data[data*convert < vmin*10.] = vmin*10./convert - im = ax_curr.imshow(convert*data, origin='lower', **kwargs) + data[data * convert < vmin * 10.0] = vmin * 10.0 / convert + im = ax_curr.imshow(convert * data, origin="lower", **kwargs) 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=1, color='grey', alpha=0.5) - ax_curr.plot([0, data.shape[1]-1], [data.shape[1]/2, data.shape[1]/2], '--', lw=1, color='grey', alpha=0.5) - ax_curr.annotate(instr+":"+rootname, color='white', fontsize=5, xy=(0.01, 1.00), xycoords='axes fraction', verticalalignment='top', horizontalalignment='left') - ax_curr.annotate(filt, color='white', fontsize=10, xy=(0.01, 0.01), xycoords='axes fraction', verticalalignment='bottom', horizontalalignment='left') - ax_curr.annotate(exptime, color='white', fontsize=5, xy=(1.00, 0.01), xycoords='axes fraction', verticalalignment='bottom', horizontalalignment='right') + ax_curr.plot([data.shape[1] / 2, data.shape[1] / 2], [0, data.shape[0] - 1], "--", lw=1, color="grey", alpha=0.5) + ax_curr.plot([0, data.shape[1] - 1], [data.shape[1] / 2, data.shape[1] / 2], "--", lw=1, color="grey", alpha=0.5) + ax_curr.annotate( + instr + ":" + rootname, color="white", fontsize=5, xy=(0.01, 1.00), xycoords="axes fraction", verticalalignment="top", horizontalalignment="left" + ) + ax_curr.annotate(filt, color="white", fontsize=10, xy=(0.01, 0.01), xycoords="axes fraction", verticalalignment="bottom", horizontalalignment="left") + ax_curr.annotate(exptime, color="white", fontsize=5, xy=(1.00, 0.01), xycoords="axes fraction", verticalalignment="bottom", horizontalalignment="right") # 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"Flux [$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"Flux [$ergs \cdot cm^{-2} \cdot s^{-1} \cdot \AA^{-1}$]") if not (savename is None): # fig.suptitle(savename) - if not savename[-4:] in ['.png', '.jpg', '.pdf']: - savename += '.pdf' - fig.savefig(path_join(plots_folder, savename), bbox_inches='tight') + if not savename[-4:] in [".png", ".jpg", ".pdf"]: + savename += ".pdf" + fig.savefig(path_join(plots_folder, savename), bbox_inches="tight") plt.show() return 0 @@ -165,10 +171,10 @@ def plot_Stokes(Stokes, savename=None, plots_folder=""): Defaults to current folder. """ # Get data - stkI = Stokes['I_stokes'].data.copy() - stkQ = Stokes['Q_stokes'].data.copy() - stkU = Stokes['U_stokes'].data.copy() - data_mask = Stokes['Data_mask'].data.astype(bool) + stkI = Stokes["I_stokes"].data.copy() + stkQ = Stokes["Q_stokes"].data.copy() + stkU = Stokes["U_stokes"].data.copy() + data_mask = Stokes["Data_mask"].data.astype(bool) for dataset in [stkI, stkQ, stkU]: dataset[np.logical_not(data_mask)] = np.nan @@ -176,36 +182,47 @@ def plot_Stokes(Stokes, savename=None, plots_folder=""): wcs = WCS(Stokes[0]).deepcopy() # Plot figure - plt.rcParams.update({'font.size': 10}) + plt.rcParams.update({"font.size": 10}) fig, (axI, axQ, axU) = plt.subplots(ncols=3, figsize=(20, 6), subplot_kw=dict(projection=wcs)) fig.subplots_adjust(hspace=0, wspace=0.75, bottom=0.01, top=0.99, left=0.08, right=0.95) fig.suptitle("I, Q, U Stokes parameters") - imI = axI.imshow(stkI, origin='lower', cmap='inferno') - fig.colorbar(imI, ax=axI, aspect=50, shrink=0.50, pad=0.025, label='counts/sec') - axI.set(xlabel="RA", ylabel='DEC', title=r"$I_{stokes}$") + imI = axI.imshow(stkI, origin="lower", cmap="inferno") + fig.colorbar(imI, ax=axI, aspect=50, shrink=0.50, pad=0.025, label="counts/sec") + axI.set(xlabel="RA", ylabel="DEC", title=r"$I_{stokes}$") - imQ = axQ.imshow(stkQ, origin='lower', cmap='inferno') - fig.colorbar(imQ, ax=axQ, aspect=50, shrink=0.50, pad=0.025, label='counts/sec') - axQ.set(xlabel="RA", ylabel='DEC', title=r"$Q_{stokes}$") + imQ = axQ.imshow(stkQ, origin="lower", cmap="inferno") + fig.colorbar(imQ, ax=axQ, aspect=50, shrink=0.50, pad=0.025, label="counts/sec") + axQ.set(xlabel="RA", ylabel="DEC", title=r"$Q_{stokes}$") - imU = axU.imshow(stkU, origin='lower', cmap='inferno') - fig.colorbar(imU, ax=axU, aspect=50, shrink=0.50, pad=0.025, label='counts/sec') - axU.set(xlabel="RA", ylabel='DEC', title=r"$U_{stokes}$") + imU = axU.imshow(stkU, origin="lower", cmap="inferno") + fig.colorbar(imU, ax=axU, aspect=50, shrink=0.50, pad=0.025, label="counts/sec") + axU.set(xlabel="RA", ylabel="DEC", title=r"$U_{stokes}$") if not (savename is None): # fig.suptitle(savename+"_IQU") - if not savename[-4:] in ['.png', '.jpg', '.pdf']: - savename += '_IQU.pdf' + if not savename[-4:] in [".png", ".jpg", ".pdf"]: + savename += "_IQU.pdf" else: - savename = savename[:-4]+"_IQU"+savename[-4:] - fig.savefig(path_join(plots_folder, savename), bbox_inches='tight') + savename = savename[:-4] + "_IQU" + savename[-4:] + fig.savefig(path_join(plots_folder, savename), bbox_inches="tight") plt.show() return 0 -def polarization_map(Stokes, data_mask=None, rectangle=None, SNRp_cut=3., SNRi_cut=3., - flux_lim=None, step_vec=1, vec_scale=2., savename=None, plots_folder="", display="default"): +def polarization_map( + Stokes, + data_mask=None, + rectangle=None, + SNRp_cut=3.0, + SNRi_cut=3.0, + flux_lim=None, + step_vec=1, + vec_scale=2.0, + savename=None, + plots_folder="", + display="default", +): """ Plots polarization map from Stokes HDUList. ---------- @@ -256,14 +273,14 @@ def polarization_map(Stokes, data_mask=None, rectangle=None, SNRp_cut=3., SNRi_c The figure and ax created for interactive contour maps. """ # Get data - stkI = Stokes['I_stokes'].data.copy() - stk_cov = Stokes['IQU_cov_matrix'].data.copy() - pol = Stokes['Pol_deg_debiased'].data.copy() - pol_err = Stokes['Pol_deg_err'].data.copy() - pang = Stokes['Pol_ang'].data.copy() + stkI = Stokes["I_stokes"].data.copy() + stk_cov = Stokes["IQU_cov_matrix"].data.copy() + pol = Stokes["Pol_deg_debiased"].data.copy() + pol_err = Stokes["Pol_deg_err"].data.copy() + pang = Stokes["Pol_ang"].data.copy() try: if data_mask is None: - data_mask = Stokes['Data_mask'].data.astype(bool).copy() + data_mask = Stokes["Data_mask"].data.astype(bool).copy() except KeyError: data_mask = np.ones(stkI.shape).astype(bool) @@ -273,23 +290,23 @@ def polarization_map(Stokes, data_mask=None, rectangle=None, SNRp_cut=3., SNRi_c for j in range(3): stk_cov[i][j][np.logical_not(data_mask)] = np.nan - pivot_wav = Stokes[0].header['photplam'] - convert_flux = Stokes[0].header['photflam'] + pivot_wav = Stokes[0].header["photplam"] + convert_flux = Stokes[0].header["photflam"] wcs = WCS(Stokes[0]).deepcopy() # Plot Stokes parameters map - if display is None or display.lower() in ['default']: + if display is None or display.lower() in ["default"]: plot_Stokes(Stokes, savename=savename, plots_folder=plots_folder) # Compute SNR and apply cuts poldata, pangdata = pol.copy(), pang.copy() maskP = pol_err > 0 - SNRp = np.ones(pol.shape)*np.nan - SNRp[maskP] = pol[maskP]/pol_err[maskP] + SNRp = np.ones(pol.shape) * np.nan + SNRp[maskP] = pol[maskP] / pol_err[maskP] maskI = stk_cov[0, 0] > 0 - SNRi = np.ones(stkI.shape)*np.nan - SNRi[maskI] = stkI[maskI]/np.sqrt(stk_cov[0, 0][maskI]) + SNRi = np.ones(stkI.shape) * np.nan + SNRi[maskI] = stkI[maskI] / np.sqrt(stk_cov[0, 0][maskI]) mask = (SNRp > SNRp_cut) * (SNRi > SNRi_cut) poldata[np.logical_not(mask)] = np.nan @@ -303,163 +320,217 @@ def polarization_map(Stokes, data_mask=None, rectangle=None, SNRp_cut=3., SNRi_c print("No pixel with polarization information above requested SNR.") # Plot the map - plt.rcParams.update({'font.size': 10}) + plt.rcParams.update({"font.size": 10}) plt.rcdefaults() - fig, ax = plt.subplots(figsize=(10, 10), layout='constrained', subplot_kw=dict(projection=wcs)) - ax.set(aspect='equal', fc='k') + fig, ax = plt.subplots(figsize=(10, 10), layout="constrained", subplot_kw=dict(projection=wcs)) + ax.set(aspect="equal", fc="k") # fig.subplots_adjust(hspace=0, wspace=0, left=0.102, right=1.02) - if display.lower() in ['intensity']: + if display.lower() in ["intensity"]: # If no display selected, show intensity map - display = 'i' + display = "i" if flux_lim is None: - if mask.sum() > 0.: - vmin, vmax = 1./2.*np.median(np.sqrt(stk_cov[0, 0][mask])*convert_flux), np.max(stkI[stkI > 0.]*convert_flux) + if mask.sum() > 0.0: + vmin, vmax = 1.0 / 2.0 * np.median(np.sqrt(stk_cov[0, 0][mask]) * convert_flux), np.max(stkI[stkI > 0.0] * convert_flux) else: - vmin, vmax = 1./2.*np.median(np.sqrt(stk_cov[0, 0][stkI > 0.])*convert_flux), np.max(stkI[stkI > 0.]*convert_flux) + 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) else: vmin, vmax = flux_lim - im = ax.imshow(stkI*convert_flux, norm=LogNorm(vmin, vmax), aspect='equal', cmap='inferno', alpha=1.) + im = ax.imshow(stkI * convert_flux, norm=LogNorm(vmin, vmax), aspect="equal", cmap="inferno", alpha=1.0) 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}$]") - levelsI = np.array([0.8, 2., 5., 10., 20., 50.])/100.*vmax + levelsI = np.array([0.8, 2.0, 5.0, 10.0, 20.0, 50.0]) / 100.0 * vmax print("Total flux contour levels : ", levelsI) - ax.contour(stkI*convert_flux, levels=levelsI, colors='grey', linewidths=0.5) - elif display.lower() in ['pol_flux']: + ax.contour(stkI * convert_flux, levels=levelsI, colors="grey", linewidths=0.5) + elif display.lower() in ["pol_flux"]: # Display polarization flux - display = 'pf' + display = "pf" if flux_lim is None: - if mask.sum() > 0.: - vmin, vmax = 1./2.*np.median(np.sqrt(stk_cov[0, 0][mask])*convert_flux), np.max(stkI[stkI > 0.]*convert_flux) + if mask.sum() > 0.0: + vmin, vmax = 1.0 / 2.0 * np.median(np.sqrt(stk_cov[0, 0][mask]) * convert_flux), np.max(stkI[stkI > 0.0] * convert_flux) else: - vmin, vmax = 1./2.*np.median(np.sqrt(stk_cov[0, 0][stkI > 0.])*convert_flux), np.max(stkI[stkI > 0.]*convert_flux) + 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) else: vmin, vmax = flux_lim - im = ax.imshow(stkI*convert_flux*pol, norm=LogNorm(vmin, vmax), aspect='equal', cmap='inferno', alpha=1.) + im = ax.imshow(stkI * convert_flux * pol, norm=LogNorm(vmin, vmax), aspect="equal", cmap="inferno", alpha=1.0) fig.colorbar(im, ax=ax, aspect=50, shrink=0.75, pad=0.025, label=r"$F_{\lambda} \cdot P$ [$ergs \cdot cm^{-2} \cdot s^{-1} \cdot \AA^{-1}$]") - levelsPf = np.linspace(vmax*0.01, vmax*0.99, 10) + levelsPf = np.linspace(vmax * 0.01, vmax * 0.99, 10) print("Polarized flux contour levels : ", levelsPf) - ax.contour(stkI*convert_flux*pol, levels=levelsPf, colors='grey', linewidths=0.5) - elif display.lower() in ['p', 'pol', 'pol_deg']: + ax.contour(stkI * convert_flux * pol, levels=levelsPf, colors="grey", linewidths=0.5) + elif display.lower() in ["p", "pol", "pol_deg"]: # Display polarization degree map - display = 'p' - vmin, vmax = 0., 100. - im = ax.imshow(pol*100., vmin=vmin, vmax=vmax, aspect='equal', cmap='inferno', alpha=1.) + display = "p" + vmin, vmax = 0.0, 100.0 + im = ax.imshow(pol * 100.0, vmin=vmin, vmax=vmax, aspect="equal", cmap="inferno", alpha=1.0) fig.colorbar(im, ax=ax, aspect=50, shrink=0.75, pad=0.025, label=r"$P$ [%]") - elif display.lower() in ['pa', 'pang', 'pol_ang']: + elif display.lower() in ["pa", "pang", "pol_ang"]: # Display polarization degree map - display = 'pa' - vmin, vmax = 0., 180. - im = ax.imshow(princ_angle(pang), vmin=vmin, vmax=vmax, aspect='equal', cmap='inferno', alpha=1.) + display = "pa" + vmin, vmax = 0.0, 180.0 + im = ax.imshow(princ_angle(pang), vmin=vmin, vmax=vmax, aspect="equal", cmap="inferno", alpha=1.0) fig.colorbar(im, ax=ax, aspect=50, shrink=0.75, pad=0.025, label=r"$\theta_P$ [°]") - elif display.lower() in ['s_p', 'pol_err', 'pol_deg_err']: + elif display.lower() in ["s_p", "pol_err", "pol_deg_err"]: # Display polarization degree error map - display = 's_p' + display = "s_p" if (SNRp > SNRp_cut).any(): - vmin, vmax = 0., np.max([pol_err[SNRp > SNRp_cut].max(), 1.])*100. - im = ax.imshow(pol_err*100., vmin=vmin, vmax=vmax, aspect='equal', cmap='inferno_r', alpha=1.) + vmin, vmax = 0.0, np.max([pol_err[SNRp > SNRp_cut].max(), 1.0]) * 100.0 + im = ax.imshow(pol_err * 100.0, vmin=vmin, vmax=vmax, aspect="equal", cmap="inferno_r", alpha=1.0) else: - vmin, vmax = 0., 100. - im = ax.imshow(pol_err*100., vmin=vmin, vmax=vmax, aspect='equal', cmap='inferno_r', alpha=1.) + 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.75, pad=0.025, label=r"$\sigma_P$ [%]") - elif display.lower() in ['s_i', 'i_err']: + elif display.lower() in ["s_i", "i_err"]: # Display intensity error map - display = 's_i' + display = "s_i" if (SNRi > SNRi_cut).any(): - vmin, vmax = 1./2.*np.median(np.sqrt(stk_cov[0, 0][stk_cov[0, 0] > 0.]) * - convert_flux), np.max(np.sqrt(stk_cov[0, 0][stk_cov[0, 0] > 0.])*convert_flux) - im = ax.imshow(np.sqrt(stk_cov[0, 0])*convert_flux, norm=LogNorm(vmin, vmax), aspect='equal', cmap='inferno_r', alpha=1.) + vmin, vmax = ( + 1.0 / 2.0 * np.median(np.sqrt(stk_cov[0, 0][stk_cov[0, 0] > 0.0]) * convert_flux), + np.max(np.sqrt(stk_cov[0, 0][stk_cov[0, 0] > 0.0]) * convert_flux), + ) + 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='inferno', alpha=1.) + im = ax.imshow(np.sqrt(stk_cov[0, 0]) * convert_flux, aspect="equal", cmap="inferno", alpha=1.0) fig.colorbar(im, ax=ax, aspect=50, shrink=0.75, pad=0.025, label=r"$\sigma_I$ [$ergs \cdot cm^{-2} \cdot s^{-1} \cdot \AA^{-1}$]") - elif display.lower() in ['snr', 'snri']: + elif display.lower() in ["snr", "snri"]: # Display I_stokes signal-to-noise map - display = 'snri' - vmin, vmax = 0., np.max(SNRi[np.isfinite(SNRi)]) - if vmax*0.99 > SNRi_cut: - im = ax.imshow(SNRi, vmin=vmin, vmax=vmax, aspect='equal', cmap='inferno', alpha=1.) - levelsSNRi = np.linspace(SNRi_cut, vmax*0.99, 5) + display = "snri" + vmin, vmax = 0.0, np.max(SNRi[np.isfinite(SNRi)]) + if vmax * 0.99 > SNRi_cut: + im = ax.imshow(SNRi, vmin=vmin, vmax=vmax, aspect="equal", cmap="inferno", alpha=1.0) + levelsSNRi = np.linspace(SNRi_cut, vmax * 0.99, 5) print("SNRi contour levels : ", levelsSNRi) - ax.contour(SNRi, levels=levelsSNRi, colors='grey', linewidths=0.5) + ax.contour(SNRi, levels=levelsSNRi, colors="grey", linewidths=0.5) else: - im = ax.imshow(SNRi, aspect='equal', cmap='inferno', alpha=1.) + im = ax.imshow(SNRi, aspect="equal", cmap="inferno", alpha=1.0) fig.colorbar(im, ax=ax, aspect=50, shrink=0.75, pad=0.025, label=r"$I_{Stokes}/\sigma_{I}$") - elif display.lower() in ['snrp']: + elif display.lower() in ["snrp"]: # Display polarization degree signal-to-noise map - display = 'snrp' - vmin, vmax = 0., np.max(SNRp[np.isfinite(SNRp)]) - if vmax*0.99 > SNRp_cut: - im = ax.imshow(SNRp, vmin=vmin, vmax=vmax, aspect='equal', cmap='inferno', alpha=1.) - levelsSNRp = np.linspace(SNRp_cut, vmax*0.99, 5) + display = "snrp" + vmin, vmax = 0.0, np.max(SNRp[np.isfinite(SNRp)]) + if vmax * 0.99 > SNRp_cut: + im = ax.imshow(SNRp, vmin=vmin, vmax=vmax, aspect="equal", cmap="inferno", alpha=1.0) + levelsSNRp = np.linspace(SNRp_cut, vmax * 0.99, 5) print("SNRp contour levels : ", levelsSNRp) - ax.contour(SNRp, levels=levelsSNRp, colors='grey', linewidths=0.5) + ax.contour(SNRp, levels=levelsSNRp, colors="grey", linewidths=0.5) else: - im = ax.imshow(SNRp, aspect='equal', cmap='inferno', alpha=1.) + im = ax.imshow(SNRp, aspect="equal", cmap="inferno", alpha=1.0) fig.colorbar(im, ax=ax, aspect=50, shrink=0.75, pad=0.025, label=r"$P/\sigma_{P}$") else: # Defaults to intensity map - if mask.sum() > 0.: - vmin, vmax = 1.*np.mean(np.sqrt(stk_cov[0, 0][mask])*convert_flux), np.max(stkI[stkI > 0.]*convert_flux) + if 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.*np.mean(np.sqrt(stk_cov[0, 0][stkI > 0.])*convert_flux), np.max(stkI[stkI > 0.]*convert_flux) - im = ax.imshow(stkI*convert_flux, norm=LogNorm(vmin, vmax), aspect='equal', cmap='inferno', alpha=1.) + 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="inferno", alpha=1.0) 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$]") # Get integrated values from header I_diluted = stkI[data_mask].sum() I_diluted_err = np.sqrt(np.sum(stk_cov[0, 0][data_mask])) - P_diluted = Stokes[0].header['P_int'] - P_diluted_err = Stokes[0].header['P_int_err'] - PA_diluted = Stokes[0].header['PA_int'] - PA_diluted_err = Stokes[0].header['PA_int_err'] + P_diluted = Stokes[0].header["P_int"] + P_diluted_err = Stokes[0].header["P_int_err"] + PA_diluted = Stokes[0].header["PA_int"] + PA_diluted_err = Stokes[0].header["PA_int_err"] - plt.rcParams.update({'font.size': 12}) - px_size = wcs.wcs.get_cdelt()[0]*3600. - px_sc = AnchoredSizeBar(ax.transData, 1./px_size, '1 arcsec', 3, pad=0.5, sep=5, borderpad=0.5, frameon=False, size_vertical=0.005, color='w') - north_dir = AnchoredDirectionArrows(ax.transAxes, "E", "N", length=-0.08, fontsize=0.025, loc=1, aspect_ratio=-1, sep_y=0.01, sep_x=0.01, back_length=0., head_length=10., head_width=10., - angle=-Stokes[0].header['orientat'], text_props={'ec': 'k', 'fc': 'w', 'alpha': 1, 'lw': -0.2}, arrow_props={'ec': 'k', 'fc': 'w', 'alpha': 1, 'lw': 1}) + plt.rcParams.update({"font.size": 12}) + 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") + north_dir = AnchoredDirectionArrows( + ax.transAxes, + "E", + "N", + length=-0.08, + fontsize=0.025, + loc=1, + aspect_ratio=-1, + sep_y=0.01, + sep_x=0.01, + back_length=0.0, + head_length=10.0, + head_width=10.0, + angle=-Stokes[0].header["orientat"], + text_props={"ec": "k", "fc": "w", "alpha": 1, "lw": -0.2}, + arrow_props={"ec": "k", "fc": "w", "alpha": 1, "lw": 1}, + ) - if display.lower() in ['i', 's_i', 'snri', 'pf', 'p', 'pa', 's_p', 'snrp']: + if display.lower() in ["i", "s_i", "snri", "pf", "p", "pa", "s_p", "snrp"]: if step_vec == 0: - poldata[np.isfinite(poldata)] = 1./2. + poldata[np.isfinite(poldata)] = 1.0 / 2.0 step_vec = 1 - vec_scale = 2. + vec_scale = 2.0 X, Y = np.meshgrid(np.arange(stkI.shape[1]), np.arange(stkI.shape[0])) - U, V = poldata*np.cos(np.pi/2.+pangdata*np.pi/180.), poldata*np.sin(np.pi/2.+pangdata*np.pi/180.) - ax.quiver(X[::step_vec, ::step_vec], Y[::step_vec, ::step_vec], U[::step_vec, ::step_vec], V[::step_vec, ::step_vec], units='xy', angles='uv', - scale=1./vec_scale, scale_units='xy', pivot='mid', headwidth=0., headlength=0., headaxislength=0., width=0.5, linewidth=0.75, color='w', edgecolor='k') - pol_sc = AnchoredSizeBar(ax.transData, vec_scale, r"$P$= 100 %", 4, pad=0.5, sep=5, borderpad=0.5, frameon=False, size_vertical=0.005, color='w') + 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], + U[::step_vec, ::step_vec], + V[::step_vec, ::step_vec], + units="xy", + angles="uv", + scale=1.0 / vec_scale, + scale_units="xy", + pivot="mid", + headwidth=0.0, + headlength=0.0, + headaxislength=0.0, + width=0.5, + linewidth=0.75, + color="w", + edgecolor="k", + ) + pol_sc = AnchoredSizeBar(ax.transData, vec_scale, 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) ax.add_artist(px_sc) 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*convert_flux, I_diluted_err*convert_flux, 2))+"\n"+r"$P^{{int}}$ = {0:.1f} $\pm$ {1:.1f} %".format(P_diluted*100., P_diluted_err * - 100.)+"\n"+r"$\theta_{{P}}^{{int}}$ = {0:.1f} $\pm$ {1:.1f} °".format(PA_diluted, PA_diluted_err), color='white', xy=(0.01, 1.00), xycoords='axes fraction', path_effects=[pe.withStroke(linewidth=0.5, foreground='k')], verticalalignment='top', horizontalalignment='left') + 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) + ) + + "\n" + + r"$P^{{int}}$ = {0:.1f} $\pm$ {1:.1f} %".format(P_diluted * 100.0, P_diluted_err * 100.0) + + "\n" + + r"$\theta_{{P}}^{{int}}$ = {0:.1f} $\pm$ {1:.1f} °".format(PA_diluted, PA_diluted_err), + color="white", + xy=(0.01, 1.00), + xycoords="axes fraction", + path_effects=[pe.withStroke(linewidth=0.5, foreground="k")], + verticalalignment="top", + horizontalalignment="left", + ) else: - if display.lower() == 'default': + if display.lower() == "default": ax.add_artist(px_sc) 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*convert_flux, I_diluted_err*convert_flux, 2)), - color='white', xy=(0.01, 1.00), xycoords='axes fraction', path_effects=[pe.withStroke(linewidth=0.5, foreground='k')], verticalalignment='top', horizontalalignment='left') + 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) + ), + color="white", + xy=(0.01, 1.00), + xycoords="axes fraction", + path_effects=[pe.withStroke(linewidth=0.5, foreground="k")], + verticalalignment="top", + horizontalalignment="left", + ) # Display instrument FOV if not (rectangle is None): x, y, width, height, angle, color = rectangle - x, y = np.array([x, y]) - np.array(stkI.shape)/2. - ax.add_patch(Rectangle((x, y), width, height, angle=angle, - edgecolor=color, fill=False)) + x, y = np.array([x, y]) - np.array(stkI.shape) / 2.0 + ax.add_patch(Rectangle((x, y), width, height, angle=angle, edgecolor=color, fill=False)) # ax.coords.grid(True, color='white', ls='dotted', alpha=0.5) - ax.coords[0].set_axislabel('Right Ascension (J2000)') - ax.coords[0].set_axislabel_position('t') - ax.coords[0].set_ticklabel_position('t') - ax.set_ylabel('Declination (J2000)', labelpad=-1) + ax.coords[0].set_axislabel("Right Ascension (J2000)") + ax.coords[0].set_axislabel_position("t") + ax.coords[0].set_ticklabel_position("t") + ax.set_ylabel("Declination (J2000)", labelpad=-1) 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=200) + if savename[-4:] not in [".png", ".jpg", ".pdf"]: + savename += ".pdf" + fig.savefig(path_join(plots_folder, savename), bbox_inches="tight", dpi=200) plt.show() return fig, ax @@ -475,8 +546,8 @@ class align_maps(object): self.map = map self.other = other_map - self.map_path = self.map.fileinfo(0)['filename'] - self.other_path = self.other.fileinfo(0)['filename'] + self.map_path = self.map.fileinfo(0)["filename"] + self.other_path = self.other.fileinfo(0)["filename"] self.map_header = fits.getheader(self.map_path) self.other_header = fits.getheader(self.other_path) @@ -495,16 +566,26 @@ class align_maps(object): elif len(self.other_data.shape) == 3: 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}$") if "PHOTFLAM" in list( - self.map_header.keys()) else (1., 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}$") if "PHOTFLAM" in list( - self.other_header.keys()) else (1., 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'] - self.other_observer = "/".join([self.other_header['telescop'], self.other_header['instrume']] - ) if "INSTRUME" in list(self.other_header.keys()) else self.other_header['telescop'] + self.map_convert, self.map_unit = ( + (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") + ) + self.other_convert, self.other_unit = ( + (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") + ) + 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"] + ) + self.other_observer = ( + "/".join([self.other_header["telescop"], self.other_header["instrume"]]) + if "INSTRUME" in list(self.other_header.keys()) + else self.other_header["telescop"] + ) - plt.rcParams.update({'font.size': 10}) + plt.rcParams.update({"font.size": 10}) fontprops = fm.FontProperties(size=16) self.fig_align = plt.figure(figsize=(20, 10)) self.map_ax = self.fig_align.add_subplot(121, projection=self.map_wcs) @@ -512,66 +593,147 @@ class align_maps(object): # Plot the UV map other_kwargs = deepcopy(kwargs) - vmin, vmax = self.map_data[self.map_data > 0.].max()/1e3*self.map_convert, self.map_data[self.map_data > 0.].max()*self.map_convert + 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: for key_i, val_i in value: kwargs[key_i] = val_i - self.map_ax.imshow(self.map_data*self.map_convert, aspect='equal', **kwargs) + self.map_ax.imshow(self.map_data * self.map_convert, aspect="equal", **kwargs) - if kwargs['cmap'] in ['inferno', 'magma', 'Greys_r', 'binary_r', 'gist_yarg_r', 'gist_gray', 'gray', 'bone', 'pink', 'hot', 'afmhot', 'gist_heat', 'copper', 'gist_earth', 'gist_stern', 'gnuplot', 'gnuplot2', 'CMRmap', 'cubehelix', 'nipy_spectral', 'gist_ncar', 'viridis']: - self.map_ax.set_facecolor('black') - self.other_ax.set_facecolor('black') + if kwargs["cmap"] in [ + "inferno", + "magma", + "Greys_r", + "binary_r", + "gist_yarg_r", + "gist_gray", + "gray", + "bone", + "pink", + "hot", + "afmhot", + "gist_heat", + "copper", + "gist_earth", + "gist_stern", + "gnuplot", + "gnuplot2", + "CMRmap", + "cubehelix", + "nipy_spectral", + "gist_ncar", + "viridis", + ]: + self.map_ax.set_facecolor("black") + self.other_ax.set_facecolor("black") font_color = "white" else: - self.map_ax.set_facecolor('white') - self.other_ax.set_facecolor('white') + self.map_ax.set_facecolor("white") + self.other_ax.set_facecolor("white") font_color = "black" - px_size1 = self.map_wcs.wcs.get_cdelt()[0]*3600. - px_sc1 = AnchoredSizeBar(self.map_ax.transData, 1./px_size1, '1 arcsec', 3, pad=0.5, sep=5, borderpad=0.5, - frameon=False, size_vertical=0.005, color=font_color, fontproperties=fontprops) + px_size1 = self.map_wcs.wcs.get_cdelt()[0] * 3600.0 + px_sc1 = AnchoredSizeBar( + self.map_ax.transData, + 1.0 / px_size1, + "1 arcsec", + 3, + pad=0.5, + sep=5, + borderpad=0.5, + frameon=False, + size_vertical=0.005, + color=font_color, + fontproperties=fontprops, + ) self.map_ax.add_artist(px_sc1) - if 'PHOTPLAM' in list(self.map_header.keys()): - self.map_ax.annotate(r"$\lambda$ = {0:.0f} $\AA$".format(self.map_header['photplam']), color=font_color, fontsize=12, xy=( - 0.01, 0.93), xycoords='axes fraction', path_effects=[pe.withStroke(linewidth=0.5, foreground='k')]) - if 'ORIENTAT' in list(self.map_header.keys()): - north_dir1 = AnchoredDirectionArrows(self.map_ax.transAxes, "E", "N", length=-0.08, fontsize=0.03, loc=1, aspect_ratio=-1, sep_y=0.01, - sep_x=0.01, angle=-self.map_header['orientat'], color=font_color, arrow_props={'ec': 'k', 'fc': 'w', 'alpha': 1, 'lw': 0.5}) + if "PHOTPLAM" in list(self.map_header.keys()): + self.map_ax.annotate( + r"$\lambda$ = {0:.0f} $\AA$".format(self.map_header["photplam"]), + color=font_color, + fontsize=12, + xy=(0.01, 0.93), + xycoords="axes fraction", + path_effects=[pe.withStroke(linewidth=0.5, foreground="k")], + ) + if "ORIENTAT" in list(self.map_header.keys()): + north_dir1 = AnchoredDirectionArrows( + self.map_ax.transAxes, + "E", + "N", + length=-0.08, + fontsize=0.03, + loc=1, + aspect_ratio=-1, + sep_y=0.01, + sep_x=0.01, + angle=-self.map_header["orientat"], + color=font_color, + arrow_props={"ec": "k", "fc": "w", "alpha": 1, "lw": 0.5}, + ) self.map_ax.add_artist(north_dir1) - self.cr_map, = self.map_ax.plot(*(self.map_wcs.wcs.crpix-(1., 1.)), 'r+') + (self.cr_map,) = self.map_ax.plot(*(self.map_wcs.wcs.crpix - (1.0, 1.0)), "r+") self.map_ax.set_title("{0:s} observation\nClick on selected point of reference.".format(self.map_observer)) self.map_ax.set_xlabel(label="Right Ascension (J2000)") self.map_ax.set_ylabel(label="Declination (J2000)", labelpad=-1) # Plot the other map - vmin, vmax = self.other_data[self.other_data > 0.].max()/1e3*self.other_convert, self.other_data[self.other_data > 0.].max()*self.other_convert + vmin, vmax = 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)]]]]: try: _ = other_kwargs[key] except KeyError: for key_i, val_i in value: other_kwargs[key_i] = val_i - self.other_ax.imshow(self.other_data*self.other_convert, aspect='equal', **other_kwargs) + self.other_ax.imshow(self.other_data * self.other_convert, aspect="equal", **other_kwargs) - px_size2 = self.other_wcs.wcs.get_cdelt()[0]*3600. - px_sc2 = AnchoredSizeBar(self.other_ax.transData, 1./px_size2, '1 arcsec', 3, pad=0.5, sep=5, borderpad=0.5, - frameon=False, size_vertical=0.005, color=font_color, fontproperties=fontprops) + px_size2 = self.other_wcs.wcs.get_cdelt()[0] * 3600.0 + px_sc2 = AnchoredSizeBar( + self.other_ax.transData, + 1.0 / px_size2, + "1 arcsec", + 3, + pad=0.5, + sep=5, + borderpad=0.5, + frameon=False, + size_vertical=0.005, + color=font_color, + fontproperties=fontprops, + ) self.other_ax.add_artist(px_sc2) - if 'PHOTPLAM' in list(self.other_header.keys()): - self.other_ax.annotate(r"$\lambda$ = {0:.0f} $\AA$".format(self.other_header['photplam']), color='white', fontsize=12, xy=( - 0.01, 0.93), xycoords='axes fraction', path_effects=[pe.withStroke(linewidth=0.5, foreground='k')]) - if 'ORIENTAT' in list(self.other_header.keys()): - north_dir2 = AnchoredDirectionArrows(self.map_ax.transAxes, "E", "N", length=-0.08, fontsize=0.03, loc=1, aspect_ratio=-1, sep_y=0.01, - sep_x=0.01, angle=-self.other_header['orientat'], color=font_color, arrow_props={'ec': 'k', 'fc': 'w', 'alpha': 1, 'lw': 0.5}) + if "PHOTPLAM" in list(self.other_header.keys()): + self.other_ax.annotate( + r"$\lambda$ = {0:.0f} $\AA$".format(self.other_header["photplam"]), + color="white", + fontsize=12, + xy=(0.01, 0.93), + xycoords="axes fraction", + path_effects=[pe.withStroke(linewidth=0.5, foreground="k")], + ) + if "ORIENTAT" in list(self.other_header.keys()): + north_dir2 = AnchoredDirectionArrows( + self.map_ax.transAxes, + "E", + "N", + length=-0.08, + fontsize=0.03, + loc=1, + aspect_ratio=-1, + sep_y=0.01, + sep_x=0.01, + angle=-self.other_header["orientat"], + color=font_color, + arrow_props={"ec": "k", "fc": "w", "alpha": 1, "lw": 0.5}, + ) self.other_ax.add_artist(north_dir2) - self.cr_other, = self.other_ax.plot(*(self.other_wcs.wcs.crpix-(1., 1.)), 'r+') + (self.cr_other,) = self.other_ax.plot(*(self.other_wcs.wcs.crpix - (1.0, 1.0)), "r+") self.other_ax.set_title("{0:s} observation\nClick on selected point of reference.".format(self.other_observer)) self.other_ax.set_xlabel(label="Right Ascension (J2000)") @@ -579,12 +741,12 @@ class align_maps(object): # Selection button self.axapply = self.fig_align.add_axes([0.80, 0.01, 0.1, 0.04]) - self.bapply = Button(self.axapply, 'Apply reference') + self.bapply = Button(self.axapply, "Apply reference") self.bapply.label.set_fontsize(8) self.axreset = self.fig_align.add_axes([0.60, 0.01, 0.1, 0.04]) - self.breset = Button(self.axreset, 'Leave as is') + self.breset = Button(self.axreset, "Leave as is") self.breset.label.set_fontsize(8) - self.enter = self.fig_align.canvas.mpl_connect('key_press_event', self.on_key) + self.enter = self.fig_align.canvas.mpl_connect("key_press_event", self.on_key) def on_key(self, event): if event.key.lower() == "enter": @@ -594,7 +756,7 @@ class align_maps(object): return self.map_wcs, self.other_wcs def onclick_ref(self, event) -> None: - if self.fig_align.canvas.manager.toolbar.mode == '': + if self.fig_align.canvas.manager.toolbar.mode == "": if (event.inaxes is not None) and (event.inaxes == self.map_ax): x = event.xdata y = event.ydata @@ -621,13 +783,13 @@ class align_maps(object): def apply_align(self, event=None): if np.array(self.cr_map.get_data()).shape == (2, 1): - self.map_wcs.wcs.crpix = np.array(self.cr_map.get_data())[:, 0]+(1., 1.) + self.map_wcs.wcs.crpix = np.array(self.cr_map.get_data())[:, 0] + (1.0, 1.0) else: - self.map_wcs.wcs.crpix = np.array(self.cr_map.get_data())+(1., 1.) + 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., 1.) + 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., 1.) + 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)) self.other_wcs.wcs.crval = self.map_wcs.wcs.crval self.fig_align.canvas.draw_idle() @@ -644,10 +806,10 @@ class align_maps(object): def align(self): self.fig_align.canvas.draw() - self.fig_align.canvas.mpl_connect('button_press_event', self.onclick_ref) + self.fig_align.canvas.mpl_connect("button_press_event", self.onclick_ref) self.bapply.on_clicked(self.apply_align) self.breset.on_clicked(self.reset_align) - self.fig_align.canvas.mpl_connect('close_event', self.on_close_align) + self.fig_align.canvas.mpl_connect("close_event", self.on_close_align) plt.show(block=True) return self.get_aligned_wcs() @@ -655,14 +817,14 @@ class align_maps(object): new_head = deepcopy(self.map_header) new_head.update(self.map_wcs.to_header()) new_hdul = fits.HDUList(fits.PrimaryHDU(self.map_data, new_head)) - new_hdul.writeto("_".join([path[:-5], suffix])+".fits", overwrite=True) + new_hdul.writeto("_".join([path[:-5], suffix]) + ".fits", overwrite=True) return 0 def write_other_to(self, path="other_map.fits", suffix="aligned", data_dir="."): new_head = deepcopy(self.other_header) new_head.update(self.other_wcs.to_header()) new_hdul = fits.HDUList(fits.PrimaryHDU(self.other_data, new_head)) - new_hdul.writeto("_".join([path[:-5], suffix])+".fits", overwrite=True) + new_hdul.writeto("_".join([path[:-5], suffix]) + ".fits", overwrite=True) return 0 def write_to(self, path1="map.fits", path2="other_map.fits", suffix="aligned", data_dir="."): @@ -677,114 +839,201 @@ 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, SNRp_cut=3., SNRi_cut=3., vec_scale=2, savename=None, **kwargs): + def overplot(self, levels=None, SNRp_cut=3.0, SNRi_cut=3.0, vec_scale=2, savename=None, **kwargs): self.Stokes_UV = self.map self.wcs_UV = self.map_wcs # Get Data - obj = self.Stokes_UV[0].header['targname'] - stkI = self.Stokes_UV['I_STOKES'].data - stk_cov = self.Stokes_UV['IQU_COV_MATRIX'].data - pol = deepcopy(self.Stokes_UV['POL_DEG_DEBIASED'].data) - pol_err = self.Stokes_UV['POL_DEG_ERR'].data - pang = self.Stokes_UV['POL_ANG'].data + obj = self.Stokes_UV[0].header["targname"] + stkI = self.Stokes_UV["I_STOKES"].data + stk_cov = self.Stokes_UV["IQU_COV_MATRIX"].data + pol = deepcopy(self.Stokes_UV["POL_DEG_DEBIASED"].data) + pol_err = self.Stokes_UV["POL_DEG_ERR"].data + pang = self.Stokes_UV["POL_ANG"].data other_data = self.other_data - self.other_convert = 1. - if self.other_unit.lower() == 'jy/beam': + self.other_convert = 1.0 + if self.other_unit.lower() == "jy/beam": self.other_unit = r"mJy/Beam" self.other_convert = 1e3 - other_freq = self.other_header['crval3'] if 'CRVAL3' in list(self.other_header.keys()) else 1. + other_freq = self.other_header["crval3"] if "CRVAL3" in list(self.other_header.keys()) else 1.0 - self.map_convert = self.Stokes_UV[0].header['photflam'] + self.map_convert = self.Stokes_UV[0].header["photflam"] # Compute SNR and apply cuts - pol[pol == 0.] = np.nan - SNRp = pol/pol_err - SNRp[np.isnan(SNRp)] = 0. + pol[pol == 0.0] = np.nan + SNRp = pol / pol_err + SNRp[np.isnan(SNRp)] = 0.0 pol[SNRp < SNRp_cut] = np.nan - SNRi = stkI/np.sqrt(stk_cov[0, 0]) - SNRi[np.isnan(SNRi)] = 0. + SNRi = stkI / np.sqrt(stk_cov[0, 0]) + SNRi[np.isnan(SNRi)] = 0.0 pol[SNRi < SNRi_cut] = np.nan - plt.rcParams.update({'font.size': 16}) + plt.rcParams.update({"font.size": 16}) self.fig_overplot, self.ax_overplot = plt.subplots(figsize=(10, 10), subplot_kw=dict(projection=self.wcs_UV)) 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 + 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: for key_i, val_i in value: kwargs[key_i] = val_i - if kwargs['cmap'] in ['inferno', 'magma', 'Greys_r', 'binary_r', 'gist_yarg_r', 'gist_gray', 'gray', 'bone', 'pink', 'hot', 'afmhot', 'gist_heat', 'copper', 'gist_earth', 'gist_stern', 'gnuplot', 'gnuplot2', 'CMRmap', 'cubehelix', 'nipy_spectral', 'gist_ncar', 'viridis']: - self.ax_overplot.set_facecolor('black') + if kwargs["cmap"] in [ + "inferno", + "magma", + "Greys_r", + "binary_r", + "gist_yarg_r", + "gist_gray", + "gray", + "bone", + "pink", + "hot", + "afmhot", + "gist_heat", + "copper", + "gist_earth", + "gist_stern", + "gnuplot", + "gnuplot2", + "CMRmap", + "cubehelix", + "nipy_spectral", + "gist_ncar", + "viridis", + ]: + self.ax_overplot.set_facecolor("black") font_color = "white" else: - self.ax_overplot.set_facecolor('white') + 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.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 = 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) + ) # Display full size polarization vectors if vec_scale is None: - self.vec_scale = 2. - pol[np.isfinite(pol)] = 1./2. + self.vec_scale = 2.0 + pol[np.isfinite(pol)] = 1.0 / 2.0 else: self.vec_scale = vec_scale 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.+pang*np.pi/180.), pol*np.sin(np.pi/2.+pang*np.pi/180.) - self.Q = self.ax_overplot.quiver(self.X[::step_vec, ::step_vec], self.Y[::step_vec, ::step_vec], self.U[::step_vec, ::step_vec], self.V[::step_vec, ::step_vec], units='xy', angles='uv', scale=1./self.vec_scale, - scale_units='xy', pivot='mid', headwidth=0., headlength=0., headaxislength=0., width=0.5, linewidth=0.75, color='white', edgecolor='black', label="{0:s} polarization map".format(self.map_observer)) + 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], + self.U[::step_vec, ::step_vec], + self.V[::step_vec, ::step_vec], + units="xy", + angles="uv", + scale=1.0 / self.vec_scale, + scale_units="xy", + pivot="mid", + headwidth=0.0, + headlength=0.0, + headaxislength=0.0, + width=0.5, + linewidth=0.75, + color="white", + edgecolor="black", + label="{0:s} polarization map".format(self.map_observer), + ) self.ax_overplot.autoscale(False) # Display other map as contours if levels is None: - levels = np.logspace(0., 1.9, 5)/100.*other_data[other_data > 0.].max() + levels = np.logspace(0.0, 1.9, 5) / 100.0 * other_data[other_data > 0.0].max() other_cont = self.ax_overplot.contour( - other_data*self.other_convert, transform=self.ax_overplot.get_transform(self.other_wcs.celestial), levels=levels*self.other_convert, colors='grey') + other_data * self.other_convert, + transform=self.ax_overplot.get_transform(self.other_wcs.celestial), + levels=levels * self.other_convert, + colors="grey", + ) self.ax_overplot.clabel(other_cont, inline=True, fontsize=5) 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 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), wrap=True) + 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 + ), + wrap=True, + ) # Display pixel scale and North direction fontprops = fm.FontProperties(size=16) - px_size = self.wcs_UV.wcs.get_cdelt()[0]*3600. - px_sc = AnchoredSizeBar(self.ax_overplot.transData, 1./px_size, '1 arcsec', 3, pad=0.5, sep=5, borderpad=0.5, - frameon=False, size_vertical=0.005, color=font_color, fontproperties=fontprops) + px_size = self.wcs_UV.wcs.get_cdelt()[0] * 3600.0 + px_sc = AnchoredSizeBar( + self.ax_overplot.transData, + 1.0 / px_size, + "1 arcsec", + 3, + pad=0.5, + sep=5, + borderpad=0.5, + frameon=False, + size_vertical=0.005, + color=font_color, + fontproperties=fontprops, + ) self.ax_overplot.add_artist(px_sc) - north_dir = AnchoredDirectionArrows(self.ax_overplot.transAxes, "E", "N", length=-0.08, fontsize=0.03, loc=1, aspect_ratio=-1, sep_y=0.01, - sep_x=0.01, angle=-self.Stokes_UV[0].header['orientat'], color=font_color, arrow_props={'ec': 'k', 'fc': 'w', 'alpha': 1, 'lw': 0.5}) + north_dir = AnchoredDirectionArrows( + self.ax_overplot.transAxes, + "E", + "N", + length=-0.08, + fontsize=0.03, + loc=1, + aspect_ratio=-1, + sep_y=0.01, + sep_x=0.01, + angle=-self.Stokes_UV[0].header["orientat"], + color=font_color, + arrow_props={"ec": "k", "fc": "w", "alpha": 1, "lw": 0.5}, + ) self.ax_overplot.add_artist(north_dir) - pol_sc = AnchoredSizeBar(self.ax_overplot.transData, self.vec_scale, r"$P$= 100%", 4, pad=0.5, sep=5, borderpad=0.5, - frameon=False, size_vertical=0.005, color=font_color, fontproperties=fontprops) + pol_sc = AnchoredSizeBar( + self.ax_overplot.transData, + self.vec_scale, + r"$P$= 100%", + 4, + pad=0.5, + sep=5, + borderpad=0.5, + frameon=False, + size_vertical=0.005, + color=font_color, + fontproperties=fontprops, + ) self.ax_overplot.add_artist(pol_sc) - self.cr_map, = self.ax_overplot.plot(*(self.map_wcs.celestial.wcs.crpix-(1., 1.)), 'r+') - self.cr_other, = self.ax_overplot.plot(*(self.other_wcs.celestial.wcs.crpix-(1., 1.)), 'g+', transform=self.ax_overplot.get_transform(self.other_wcs)) + (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) + ) 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) + 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 + ) labels.append("{0:s} contour".format(self.other_observer)) handles.append(Rectangle((0, 0), 1, 1, fill=False, lw=2, ec=other_cont.collections[0].get_edgecolor()[0])) - self.legend = self.ax_overplot.legend(handles=handles, labels=labels, bbox_to_anchor=( - 0., 1.02, 1., .102), loc='lower left', mode="expand", borderaxespad=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 + ) if not (savename is None): - if not savename[-4:] in ['.png', '.jpg', '.pdf']: - savename += '.pdf' - self.fig_overplot.savefig(savename, bbox_inches='tight', dpi=200) + if not savename[-4:] in [".png", ".jpg", ".pdf"]: + savename += ".pdf" + self.fig_overplot.savefig(savename, bbox_inches="tight", dpi=200) self.fig_overplot.canvas.draw() - def plot(self, levels=None, SNRp_cut=3., SNRi_cut=3., savename=None, **kwargs) -> None: + def plot(self, levels=None, SNRp_cut=3.0, SNRi_cut=3.0, savename=None, **kwargs) -> None: while not self.aligned: self.align() self.overplot(levels=levels, SNRp_cut=SNRp_cut, SNRi_cut=SNRi_cut, savename=savename, **kwargs) @@ -797,16 +1046,16 @@ 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, SNRp_cut=3., SNRi_cut=3., vec_scale=2, zoom=1, savename=None, **kwargs): + def overplot(self, levels=None, SNRp_cut=3.0, SNRi_cut=3.0, vec_scale=2, zoom=1, savename=None, **kwargs): self.Stokes_UV = self.map self.wcs_UV = self.map_wcs # Get Data - obj = self.Stokes_UV[0].header['targname'] - stkI = self.Stokes_UV['I_STOKES'].data - stk_cov = self.Stokes_UV['IQU_COV_MATRIX'].data - pol = deepcopy(self.Stokes_UV['POL_DEG_DEBIASED'].data) - pol_err = self.Stokes_UV['POL_DEG_ERR'].data - pang = self.Stokes_UV['POL_ANG'].data + obj = self.Stokes_UV[0].header["targname"] + stkI = self.Stokes_UV["I_STOKES"].data + stk_cov = self.Stokes_UV["IQU_COV_MATRIX"].data + pol = deepcopy(self.Stokes_UV["POL_DEG_DEBIASED"].data) + pol_err = self.Stokes_UV["POL_DEG_ERR"].data + pang = self.Stokes_UV["POL_ANG"].data other_data = deepcopy(self.other_data) other_wcs = self.other_wcs.deepcopy() @@ -814,96 +1063,176 @@ class overplot_chandra(align_maps): other_data = sc_zoom(other_data, zoom) other_wcs.wcs.crpix *= zoom other_wcs.wcs.cdelt /= zoom - self.other_unit = 'counts' + self.other_unit = "counts" # Compute SNR and apply cuts - pol[pol == 0.] = np.nan - SNRp = pol/pol_err - SNRp[np.isnan(SNRp)] = 0. + pol[pol == 0.0] = np.nan + SNRp = pol / pol_err + SNRp[np.isnan(SNRp)] = 0.0 pol[SNRp < SNRp_cut] = np.nan - SNRi = stkI/np.sqrt(stk_cov[0, 0]) - SNRi[np.isnan(SNRi)] = 0. + SNRi = stkI / np.sqrt(stk_cov[0, 0]) + SNRi[np.isnan(SNRi)] = 0.0 pol[SNRi < SNRi_cut] = np.nan - plt.rcParams.update({'font.size': 16}) + plt.rcParams.update({"font.size": 16}) self.fig_overplot, self.ax_overplot = plt.subplots(figsize=(11, 10), subplot_kw=dict(projection=self.wcs_UV)) 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 + 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: for key_i, val_i in value: kwargs[key_i] = val_i - if kwargs['cmap'] in ['inferno', 'magma', 'Greys_r', 'binary_r', 'gist_yarg_r', 'gist_gray', 'gray', 'bone', 'pink', 'hot', 'afmhot', 'gist_heat', 'copper', 'gist_earth', 'gist_stern', 'gnuplot', 'gnuplot2', 'CMRmap', 'cubehelix', 'nipy_spectral', 'gist_ncar', 'viridis']: - self.ax_overplot.set_facecolor('black') + if kwargs["cmap"] in [ + "inferno", + "magma", + "Greys_r", + "binary_r", + "gist_yarg_r", + "gist_gray", + "gray", + "bone", + "pink", + "hot", + "afmhot", + "gist_heat", + "copper", + "gist_earth", + "gist_stern", + "gnuplot", + "gnuplot2", + "CMRmap", + "cubehelix", + "nipy_spectral", + "gist_ncar", + "viridis", + ]: + self.ax_overplot.set_facecolor("black") font_color = "white" else: - self.ax_overplot.set_facecolor('white') + self.ax_overplot.set_facecolor("white") 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 = 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) + ) # Display full size polarization vectors if vec_scale is None: - self.vec_scale = 2. - pol[np.isfinite(pol)] = 1./2. + self.vec_scale = 2.0 + pol[np.isfinite(pol)] = 1.0 / 2.0 else: self.vec_scale = vec_scale 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.+pang*np.pi/180.), pol*np.sin(np.pi/2.+pang*np.pi/180.) - self.Q = self.ax_overplot.quiver(self.X[::step_vec, ::step_vec], self.Y[::step_vec, ::step_vec], self.U[::step_vec, ::step_vec], self.V[::step_vec, ::step_vec], units='xy', angles='uv', scale=1./self.vec_scale, - scale_units='xy', pivot='mid', headwidth=0., headlength=0., headaxislength=0., width=0.5, linewidth=0.75, color='white', edgecolor='black', label="{0:s} polarization map".format(self.map_observer)) + 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], + self.U[::step_vec, ::step_vec], + self.V[::step_vec, ::step_vec], + units="xy", + angles="uv", + scale=1.0 / self.vec_scale, + scale_units="xy", + pivot="mid", + headwidth=0.0, + headlength=0.0, + headaxislength=0.0, + width=0.5, + linewidth=0.75, + color="white", + edgecolor="black", + label="{0:s} polarization map".format(self.map_observer), + ) self.ax_overplot.autoscale(False) # Display other map as contours if levels is None: - levels = np.logspace(np.log(3)/np.log(10), 2., 5)/100.*other_data[other_data > 0.].max()*self.other_convert + levels = np.logspace(np.log(3) / np.log(10), 2.0, 5) / 100.0 * other_data[other_data > 0.0].max() * self.other_convert 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') + 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" + ) 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) + 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 + ) # Display pixel scale and North direction fontprops = fm.FontProperties(size=16) - px_size = self.wcs_UV.wcs.get_cdelt()[0]*3600. - px_sc = AnchoredSizeBar(self.ax_overplot.transData, 1./px_size, '1 arcsec', 3, pad=0.5, sep=5, borderpad=0.5, - frameon=False, size_vertical=0.005, color=font_color, fontproperties=fontprops) + px_size = self.wcs_UV.wcs.get_cdelt()[0] * 3600.0 + px_sc = AnchoredSizeBar( + self.ax_overplot.transData, + 1.0 / px_size, + "1 arcsec", + 3, + pad=0.5, + sep=5, + borderpad=0.5, + frameon=False, + size_vertical=0.005, + color=font_color, + fontproperties=fontprops, + ) self.ax_overplot.add_artist(px_sc) - north_dir = AnchoredDirectionArrows(self.ax_overplot.transAxes, "E", "N", length=-0.08, fontsize=0.03, loc=1, aspect_ratio=-1, sep_y=0.01, - sep_x=0.01, angle=-self.Stokes_UV[0].header['orientat'], color=font_color, arrow_props={'ec': 'k', 'fc': 'w', 'alpha': 1, 'lw': 0.5}) + north_dir = AnchoredDirectionArrows( + self.ax_overplot.transAxes, + "E", + "N", + length=-0.08, + fontsize=0.03, + loc=1, + aspect_ratio=-1, + sep_y=0.01, + sep_x=0.01, + angle=-self.Stokes_UV[0].header["orientat"], + color=font_color, + arrow_props={"ec": "k", "fc": "w", "alpha": 1, "lw": 0.5}, + ) self.ax_overplot.add_artist(north_dir) - pol_sc = AnchoredSizeBar(self.ax_overplot.transData, self.vec_scale, r"$P$= 100%", 4, pad=0.5, sep=5, borderpad=0.5, - frameon=False, size_vertical=0.005, color=font_color, fontproperties=fontprops) + pol_sc = AnchoredSizeBar( + self.ax_overplot.transData, + self.vec_scale, + r"$P$= 100%", + 4, + pad=0.5, + sep=5, + borderpad=0.5, + frameon=False, + size_vertical=0.005, + color=font_color, + fontproperties=fontprops, + ) self.ax_overplot.add_artist(pol_sc) - self.cr_map, = self.ax_overplot.plot(*(self.map_wcs.celestial.wcs.crpix-(1., 1.)), 'r+') - self.cr_other, = self.ax_overplot.plot(*(other_wcs.celestial.wcs.crpix-(1., 1.)), 'g+', transform=self.ax_overplot.get_transform(other_wcs)) + (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)) 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) + 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 + ) 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.collections[0].get_edgecolor()[0])) - self.legend = self.ax_overplot.legend(handles=handles, labels=labels, bbox_to_anchor=( - 0., 1.02, 1., .102), loc='lower left', mode="expand", borderaxespad=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 + ) if not (savename is None): - if not savename[-4:] in ['.png', '.jpg', '.pdf']: - savename += '.pdf' - self.fig_overplot.savefig(savename, bbox_inches='tight', dpi=200) + if not savename[-4:] in [".png", ".jpg", ".pdf"]: + savename += ".pdf" + self.fig_overplot.savefig(savename, bbox_inches="tight", dpi=200) self.fig_overplot.canvas.draw() - def plot(self, levels=None, SNRp_cut=3., SNRi_cut=3., zoom=1, savename=None, **kwargs) -> None: + def plot(self, levels=None, SNRp_cut=3.0, SNRi_cut=3.0, zoom=1, savename=None, **kwargs) -> None: while not self.aligned: self.align() self.overplot(levels=levels, SNRp_cut=SNRp_cut, SNRi_cut=SNRi_cut, zoom=zoom, savename=savename, **kwargs) @@ -916,143 +1245,229 @@ class overplot_pol(align_maps): Inherit from class align_maps in order to get the same WCS on both maps. """ - def overplot(self, levels=None, SNRp_cut=3., SNRi_cut=3., vec_scale=2., savename=None, **kwargs): + def overplot(self, levels=None, SNRp_cut=3.0, SNRi_cut=3.0, vec_scale=2.0, savename=None, **kwargs): self.Stokes_UV = self.map self.wcs_UV = self.map_wcs # Get Data - obj = self.Stokes_UV[0].header['targname'] - stkI = self.Stokes_UV['I_STOKES'].data - stk_cov = self.Stokes_UV['IQU_COV_MATRIX'].data - pol = deepcopy(self.Stokes_UV['POL_DEG_DEBIASED'].data) - pol_err = self.Stokes_UV['POL_DEG_ERR'].data - pang = self.Stokes_UV['POL_ANG'].data + obj = self.Stokes_UV[0].header["targname"] + stkI = self.Stokes_UV["I_STOKES"].data + stk_cov = self.Stokes_UV["IQU_COV_MATRIX"].data + pol = deepcopy(self.Stokes_UV["POL_DEG_DEBIASED"].data) + pol_err = self.Stokes_UV["POL_DEG_ERR"].data + pang = self.Stokes_UV["POL_ANG"].data other_data = self.other_data # Compute SNR and apply cuts - pol[pol == 0.] = np.nan - SNRp = pol/pol_err - SNRp[np.isnan(SNRp)] = 0. + pol[pol == 0.0] = np.nan + SNRp = pol / pol_err + SNRp[np.isnan(SNRp)] = 0.0 pol[SNRp < SNRp_cut] = np.nan - SNRi = stkI/np.sqrt(stk_cov[0, 0]) - SNRi[np.isnan(SNRi)] = 0. + SNRi = stkI / np.sqrt(stk_cov[0, 0]) + SNRi[np.isnan(SNRi)] = 0.0 pol[SNRi < SNRi_cut] = np.nan - plt.rcParams.update({'font.size': 16}) + plt.rcParams.update({"font.size": 16}) self.fig_overplot, self.ax_overplot = plt.subplots(figsize=(11, 10), 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) - self.fig_overplot.suptitle("{0:s} observation from {1:s} overplotted with polarization vectors and Stokes I contours from {2:s}".format( - obj, self.other_observer, self.map_observer), wrap=True) + self.fig_overplot.suptitle( + "{0:s} observation from {1:s} overplotted with polarization vectors and Stokes I contours from {2:s}".format( + obj, self.other_observer, self.map_observer + ), + wrap=True, + ) # Display "other" intensity map - vmin, vmax = other_data[other_data > 0.].max()/1e3*self.other_convert, other_data[other_data > 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]]]]: try: _ = kwargs[key] except KeyError: for key_i, val_i in value: kwargs[key_i] = val_i - if kwargs['cmap'] in ['inferno', 'magma', 'Greys_r', 'binary_r', 'gist_yarg_r', 'gist_gray', 'gray', 'bone', 'pink', 'hot', 'afmhot', 'gist_heat', 'copper', 'gist_earth', 'gist_stern', 'gnuplot', 'gnuplot2', 'CMRmap', 'cubehelix', 'nipy_spectral', 'gist_ncar', 'viridis']: - self.ax_overplot.set_facecolor('black') + if kwargs["cmap"] in [ + "inferno", + "magma", + "Greys_r", + "binary_r", + "gist_yarg_r", + "gist_gray", + "gray", + "bone", + "pink", + "hot", + "afmhot", + "gist_heat", + "copper", + "gist_earth", + "gist_stern", + "gnuplot", + "gnuplot2", + "CMRmap", + "cubehelix", + "nipy_spectral", + "gist_ncar", + "viridis", + ]: + self.ax_overplot.set_facecolor("black") font_color = "white" else: - self.ax_overplot.set_facecolor('white') + self.ax_overplot.set_facecolor("white") font_color = "black" - self.im = self.ax_overplot.imshow(other_data*self.other_convert, alpha=1., label="{0:s} observation".format(self.other_observer), **kwargs) - 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 = self.ax_overplot.imshow(other_data * self.other_convert, alpha=1.0, label="{0:s} observation".format(self.other_observer), **kwargs) + 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) + ) # Display full size polarization vectors if vec_scale is None: - self.vec_scale = 2. - pol[np.isfinite(pol)] = 1./2. + self.vec_scale = 2.0 + pol[np.isfinite(pol)] = 1.0 / 2.0 else: self.vec_scale = vec_scale step_vec = 1 - px_scale = np.abs(self.wcs_UV.wcs.get_cdelt()[0]/self.other_wcs.wcs.get_cdelt()[0]) + px_scale = np.abs(self.wcs_UV.wcs.get_cdelt()[0] / self.other_wcs.wcs.get_cdelt()[0]) 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.+pang*np.pi/180.), pol*np.sin(np.pi/2.+pang*np.pi/180.) - self.Q = self.ax_overplot.quiver(self.X[::step_vec, ::step_vec], self.Y[::step_vec, ::step_vec], self.U[::step_vec, ::step_vec], self.V[::step_vec, ::step_vec], units='xy', angles='uv', scale=px_scale/self.vec_scale, scale_units='xy', pivot='mid', - headwidth=0., headlength=0., headaxislength=0., width=0.5, linewidth=0.75, color='white', edgecolor='black', transform=self.ax_overplot.get_transform(self.wcs_UV), label="{0:s} polarization map".format(self.map_observer)) + 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], + self.U[::step_vec, ::step_vec], + self.V[::step_vec, ::step_vec], + units="xy", + angles="uv", + scale=px_scale / self.vec_scale, + scale_units="xy", + pivot="mid", + headwidth=0.0, + headlength=0.0, + headaxislength=0.0, + width=0.5, + linewidth=0.75, + color="white", + edgecolor="black", + transform=self.ax_overplot.get_transform(self.wcs_UV), + label="{0:s} polarization map".format(self.map_observer), + ) # Display Stokes I as contours if levels is None: - levels = np.array([2., 5., 10., 20., 90.])/100.*np.max(stkI[stkI > 0.])*self.map_convert - cont_stkI = self.ax_overplot.contour(stkI*self.map_convert, levels=levels, colors='grey', alpha=0.75, - transform=self.ax_overplot.get_transform(self.wcs_UV)) + 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_stkI = self.ax_overplot.contour( + 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_stkI, inline=True, fontsize=5) # Display pixel scale and North direction fontprops = fm.FontProperties(size=16) - px_size = self.other_wcs.wcs.get_cdelt()[0]*3600. - px_sc = AnchoredSizeBar(self.ax_overplot.transData, 1./px_size, '1 arcsec', 3, pad=0.5, sep=5, borderpad=0.5, - frameon=False, size_vertical=0.005, color=font_color, fontproperties=fontprops) + px_size = self.other_wcs.wcs.get_cdelt()[0] * 3600.0 + px_sc = AnchoredSizeBar( + self.ax_overplot.transData, + 1.0 / px_size, + "1 arcsec", + 3, + pad=0.5, + sep=5, + borderpad=0.5, + frameon=False, + size_vertical=0.005, + color=font_color, + fontproperties=fontprops, + ) self.ax_overplot.add_artist(px_sc) - north_dir = AnchoredDirectionArrows(self.ax_overplot.transAxes, "E", "N", length=-0.08, fontsize=0.03, loc=1, aspect_ratio=-1, sep_y=0.01, - sep_x=0.01, angle=-self.Stokes_UV[0].header['orientat'], color=font_color, arrow_props={'ec': 'k', 'fc': 'w', 'alpha': 1, 'lw': 0.5}) + north_dir = AnchoredDirectionArrows( + self.ax_overplot.transAxes, + "E", + "N", + length=-0.08, + fontsize=0.03, + loc=1, + aspect_ratio=-1, + sep_y=0.01, + sep_x=0.01, + angle=-self.Stokes_UV[0].header["orientat"], + color=font_color, + arrow_props={"ec": "k", "fc": "w", "alpha": 1, "lw": 0.5}, + ) self.ax_overplot.add_artist(north_dir) - pol_sc = AnchoredSizeBar(self.ax_overplot.transData, self.vec_scale/px_scale, r"$P$= 100%", 4, pad=0.5, sep=5, - borderpad=0.5, frameon=False, size_vertical=0.005, color=font_color, fontproperties=fontprops) + pol_sc = AnchoredSizeBar( + self.ax_overplot.transData, + self.vec_scale / px_scale, + r"$P$= 100%", + 4, + pad=0.5, + sep=5, + borderpad=0.5, + frameon=False, + size_vertical=0.005, + color=font_color, + fontproperties=fontprops, + ) self.ax_overplot.add_artist(pol_sc) - self.cr_map, = self.ax_overplot.plot(*(self.map_wcs.celestial.wcs.crpix-(1., 1.)), 'r+', transform=self.ax_overplot.get_transform(self.wcs_UV)) - self.cr_other, = self.ax_overplot.plot(*(self.other_wcs.celestial.wcs.crpix-(1., 1.)), 'g+') + (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()): - self.legend_title = r"{0:s} image at $\lambda$ = {1:.0f} $\AA$".format(self.other_observer, float(self.other_header['photplam'])) + self.legend_title = r"{0:s} image at $\lambda$ = {1:.0f} $\AA$".format(self.other_observer, float(self.other_header["photplam"])) elif "CRVAL3" in list(self.other_header.keys()): - self.legend_title = "{0:s} image at {1:.2f} GHz".format(self.other_observer, float(self.other_header['crval3'])*1e-9) + self.legend_title = "{0:s} image at {1:.2f} GHz".format(self.other_observer, float(self.other_header["crval3"]) * 1e-9) else: self.legend_title = r"{0:s} image".format(self.other_observer) 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) + 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 + ) labels.append("{0:s} Stokes I contour".format(self.map_observer)) handles.append(Rectangle((0, 0), 1, 1, fill=False, ec=cont_stkI.collections[0].get_edgecolor()[0])) - self.legend = self.ax_overplot.legend(handles=handles, labels=labels, bbox_to_anchor=( - 0., 1.02, 1., .102), loc='lower left', mode="expand", borderaxespad=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 + ) if not (savename is None): - if not savename[-4:] in ['.png', '.jpg', '.pdf']: - savename += '.pdf' - self.fig_overplot.savefig(savename, bbox_inches='tight', dpi=200) + if not savename[-4:] in [".png", ".jpg", ".pdf"]: + savename += ".pdf" + self.fig_overplot.savefig(savename, bbox_inches="tight", dpi=200) self.fig_overplot.canvas.draw() - def plot(self, levels=None, SNRp_cut=3., SNRi_cut=3., vec_scale=2., savename=None, **kwargs) -> None: + def plot(self, levels=None, SNRp_cut=3.0, SNRi_cut=3.0, vec_scale=2.0, savename=None, **kwargs) -> None: while not self.aligned: self.align() self.overplot(levels=levels, SNRp_cut=SNRp_cut, SNRi_cut=SNRi_cut, vec_scale=vec_scale, savename=savename, **kwargs) plt.show(block=True) - def add_vector(self, position='center', pol_deg=1., pol_ang=0., **kwargs): - if position == 'center': - position = np.array(self.X.shape)/2. + def add_vector(self, position="center", pol_deg=1.0, pol_ang=0.0, **kwargs): + if position == "center": + position = np.array(self.X.shape) / 2.0 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.), pol_deg*np.sin(np.radians(pol_ang)+np.pi/2.) - for key, value in [["scale", [["scale", self.vec_scale]]], ["width", [["width", 0.1]]], ["color", [["color", 'k']]]]: + 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", self.vec_scale]]], ["width", [["width", 0.1]]], ["color", [["color", "k"]]]]: try: _ = kwargs[key] except KeyError: for key_i, val_i in value: kwargs[key_i] = val_i - new_vec = self.ax_overplot.quiver(*position, u, v, units='xy', angles='uv', scale_units='xy', - pivot='mid', headwidth=0., headlength=0., headaxislength=0., **kwargs) + 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 + ) self.legend.remove() - self.legend = self.ax_overplot.legend(title=self.legend_title, bbox_to_anchor=(0., 1.02, 1., .102), loc='lower left', mode="expand", borderaxespad=0.) + self.legend = self.ax_overplot.legend( + title=self.legend_title, 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 class align_pol(object): def __init__(self, maps, **kwargs): - order = np.argsort(np.array([curr[0].header['mjd-obs'] for curr in maps])) + order = np.argsort(np.array([curr[0].header["mjd-obs"] for curr in maps])) maps = np.array(maps)[order] self.ref_map, self.other_maps = maps[0], maps[1:] @@ -1063,39 +1478,43 @@ class align_pol(object): self.kwargs = kwargs - def single_plot(self, curr_map, wcs, v_lim=None, ax_lim=None, SNRp_cut=3., SNRi_cut=3., savename=None, **kwargs): + def single_plot(self, curr_map, wcs, v_lim=None, ax_lim=None, SNRp_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 - pol = deepcopy(curr_map['POL_DEG_DEBIASED'].data) - pol_err = curr_map['POL_DEG_ERR'].data - pang = curr_map['POL_ANG'].data + stkI = curr_map["I_STOKES"].data + stk_cov = curr_map["IQU_COV_MATRIX"].data + pol = deepcopy(curr_map["POL_DEG_DEBIASED"].data) + pol_err = curr_map["POL_DEG_ERR"].data + pang = curr_map["POL_ANG"].data try: - data_mask = curr_map['DATA_MASK'].data.astype(bool) + data_mask = curr_map["DATA_MASK"].data.astype(bool) except KeyError: data_mask = np.ones(stkI.shape).astype(bool) - convert_flux = curr_map[0].header['photflam'] + convert_flux = curr_map[0].header["photflam"] # Compute SNR and apply cuts - maskpol = np.logical_and(pol_err > 0., data_mask) + maskpol = np.logical_and(pol_err > 0.0, data_mask) SNRp = np.zeros(pol.shape) - SNRp[maskpol] = pol[maskpol]/pol_err[maskpol] + SNRp[maskpol] = pol[maskpol] / pol_err[maskpol] maskI = np.logical_and(stk_cov[0, 0] > 0, data_mask) SNRi = np.zeros(stkI.shape) - SNRi[maskI] = stkI[maskI]/np.sqrt(stk_cov[0, 0][maskI]) + SNRi[maskI] = stkI[maskI] / np.sqrt(stk_cov[0, 0][maskI]) - mask = (SNRp > SNRp_cut) * (SNRi > SNRi_cut) * (pol >= 0.) + mask = (SNRp > SNRp_cut) * (SNRi > SNRi_cut) * (pol >= 0.0) pol[mask] = np.nan # Plot the map - plt.rcParams.update({'font.size': 10}) + plt.rcParams.update({"font.size": 10}) plt.rcdefaults() fig = plt.figure(figsize=(10, 10)) ax = fig.add_subplot(111, projection=wcs) - ax.set(xlabel="Right Ascension (J2000)", ylabel="Declination (J2000)", facecolor='k', - title="target {0:s} observed on {1:s}".format(curr_map[0].header['targname'], curr_map[0].header['date-obs'])) + ax.set( + xlabel="Right Ascension (J2000)", + ylabel="Declination (J2000)", + facecolor="k", + title="target {0:s} observed on {1:s}".format(curr_map[0].header["targname"], curr_map[0].header["date-obs"]), + ) fig.subplots_adjust(hspace=0, wspace=0, right=0.102) if ax_lim is not None: @@ -1104,9 +1523,9 @@ class align_pol(object): ax.set(xlim=x_lim, ylim=y_lim) if v_lim is None: - vmin, vmax = 0., np.max(stkI[stkI > 0.]*convert_flux) + vmin, vmax = 0.0, np.max(stkI[stkI > 0.0] * convert_flux) else: - vmin, vmax = v_lim*convert_flux + vmin, vmax = v_lim * convert_flux for key, value in [["cmap", [["cmap", "inferno"]]], ["norm", [["vmin", vmin], ["vmax", vmax]]]]: try: @@ -1117,33 +1536,70 @@ class align_pol(object): for key_i, val_i in value: kwargs[key_i] = val_i - im = ax.imshow(stkI*convert_flux, aspect='equal', **kwargs) + 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}$]") - px_size = wcs.wcs.get_cdelt()[0]*3600. - px_sc = AnchoredSizeBar(ax.transData, 1./px_size, '1 arcsec', 3, pad=0.5, sep=5, borderpad=0.5, frameon=False, size_vertical=0.005, color='w') + 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") ax.add_artist(px_sc) - north_dir = AnchoredDirectionArrows(ax.transAxes, "E", "N", length=-0.08, fontsize=0.025, loc=1, aspect_ratio=-1, sep_y=0.01, sep_x=0.01, back_length=0., head_length=10., head_width=10., - angle=curr_map[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}) + north_dir = AnchoredDirectionArrows( + ax.transAxes, + "E", + "N", + length=-0.08, + fontsize=0.025, + loc=1, + aspect_ratio=-1, + sep_y=0.01, + sep_x=0.01, + back_length=0.0, + head_length=10.0, + head_width=10.0, + angle=curr_map[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}, + ) ax.add_artist(north_dir) 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.+pang*np.pi/180.), pol*np.sin(np.pi/2.+pang*np.pi/180.) - ax.quiver(X[::step_vec, ::step_vec], Y[::step_vec, ::step_vec], U[::step_vec, ::step_vec], V[::step_vec, ::step_vec], units='xy', - angles='uv', scale=0.5, scale_units='xy', pivot='mid', headwidth=0., headlength=0., headaxislength=0., width=0.5, linewidth=0.75, color='w') - pol_sc = AnchoredSizeBar(ax.transData, 2., r"$P$= 100 %", 4, pad=0.5, sep=5, borderpad=0.5, frameon=False, size_vertical=0.005, color='w') + 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], + U[::step_vec, ::step_vec], + V[::step_vec, ::step_vec], + units="xy", + angles="uv", + scale=0.5, + scale_units="xy", + pivot="mid", + headwidth=0.0, + headlength=0.0, + headaxislength=0.0, + width=0.5, + 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") ax.add_artist(pol_sc) - if 'PHOTPLAM' in list(curr_map[0].header.keys()): - ax.annotate(r"$\lambda$ = {0:.0f} $\AA$".format(curr_map[0].header['photplam']), color='white', fontsize=12, xy=( - 0.01, 0.93), xycoords='axes fraction', path_effects=[pe.withStroke(linewidth=0.5, foreground='k')]) + if "PHOTPLAM" in list(curr_map[0].header.keys()): + ax.annotate( + r"$\lambda$ = {0:.0f} $\AA$".format(curr_map[0].header["photplam"]), + color="white", + fontsize=12, + xy=(0.01, 0.93), + xycoords="axes fraction", + path_effects=[pe.withStroke(linewidth=0.5, foreground="k")], + ) if savename is not None: - if savename[-4:] not in ['.png', '.jpg', '.pdf']: - savename += '.pdf' - fig.savefig(savename, bbox_inches='tight', dpi=300) + if savename[-4:] not in [".png", ".jpg", ".pdf"]: + savename += ".pdf" + fig.savefig(savename, bbox_inches="tight", dpi=300) plt.show(block=True) return fig, ax @@ -1154,27 +1610,60 @@ class align_pol(object): self.wcs, self.wcs_other[i] = curr_align.align() self.aligned[i] = curr_align.aligned - def plot(self, SNRp_cut=3., SNRi_cut=3., savename=None, **kwargs): + def plot(self, SNRp_cut=3.0, SNRi_cut=3.0, savename=None, **kwargs): while not self.aligned.all(): self.align() eps = 1e-35 - vmin = 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)]) for curr_map in self.other_maps])/2.5 - 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)]) for curr_map in self.other_maps]) - vmin = np.min([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)])])/2.5 - vmax = np.max([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)])]) + vmin = ( + 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)] + ) + for curr_map in self.other_maps + ] + ) + / 2.5 + ) + 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)]) + for curr_map in self.other_maps + ] + ) + vmin = ( + np.min( + [ + 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) + ] + ), + ] + ) + / 2.5 + ) + vmax = np.max( + [ + 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) + ] + ), + ] + ) v_lim = np.array([vmin, vmax]) - fig, ax = self.single_plot(self.ref_map, self.wcs, v_lim=v_lim, SNRp_cut=SNRp_cut, SNRi_cut=SNRi_cut, savename=savename+'_0', **kwargs) + fig, ax = self.single_plot(self.ref_map, self.wcs, v_lim=v_lim, SNRp_cut=SNRp_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, SNRp_cut=SNRp_cut, - SNRi_cut=SNRi_cut, savename=savename+'_'+str(i+1), **kwargs) + self.single_plot( + curr_map, self.wcs_other[i], v_lim=v_lim, ax_lim=ax_lim, SNRp_cut=SNRp_cut, SNRi_cut=SNRi_cut, savename=savename + "_" + str(i + 1), **kwargs + ) class crop_map(object): @@ -1191,16 +1680,16 @@ class crop_map(object): self.data = deepcopy(self.hdul[0].data) try: - self.map_convert = self.header['photflam'] + self.map_convert = self.header["photflam"] except KeyError: - self.map_convert = 1. + self.map_convert = 1.0 try: self.kwargs = kwargs except AttributeError: self.kwargs = {} # Plot the map - plt.rcParams.update({'font.size': 12}) + plt.rcParams.update({"font.size": 12}) if fig is None: self.fig = plt.figure(figsize=(15, 15)) self.fig.suptitle("Click and drag to crop to desired Region of Interest.") @@ -1208,23 +1697,22 @@ class crop_map(object): self.fig = fig if ax is None: self.ax = self.fig.add_subplot(111, projection=self.wcs) - self.mask_alpha = 1. + self.mask_alpha = 1.0 # Selection button self.axapply = self.fig.add_axes([0.80, 0.01, 0.1, 0.04]) - self.bapply = Button(self.axapply, 'Apply') + self.bapply = Button(self.axapply, "Apply") self.axreset = self.fig.add_axes([0.60, 0.01, 0.1, 0.04]) - self.breset = Button(self.axreset, 'Reset') + self.breset = Button(self.axreset, "Reset") self.embedded = False else: self.ax = ax self.mask_alpha = 0.75 - self.rect_selector = RectangleSelector(self.ax, self.onselect_crop, - button=[1]) + self.rect_selector = RectangleSelector(self.ax, self.onselect_crop, button=[1]) self.embedded = True self.display(self.data, self.wcs, self.map_convert, **self.kwargs) - self.extent = np.array([0., self.data.shape[0], 0., self.data.shape[1]]) - self.center = np.array(self.data.shape)/2 + self.extent = np.array([0.0, self.data.shape[0], 0.0, self.data.shape[1]]) + self.center = np.array(self.data.shape) / 2 self.RSextent = deepcopy(self.extent) self.RScenter = deepcopy(self.center) @@ -1240,20 +1728,26 @@ class crop_map(object): else: kwargs = {**self.kwargs, **kwargs} - vmin, vmax = np.min(data[data > 0.]*convert_flux), np.max(data[data > 0.]*convert_flux) - for key, value in [["cmap", [["cmap", "inferno"]]], ["origin", [["origin", "lower"]]], ["aspect", [["aspect", "equal"]]], ["alpha", [["alpha", self.mask_alpha]]], ["norm", [["vmin", vmin], ["vmax", vmax]]]]: + 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"]]], + ["aspect", [["aspect", "equal"]]], + ["alpha", [["alpha", self.mask_alpha]]], + ["norm", [["vmin", vmin], ["vmax", vmax]]], + ]: try: _ = kwargs[key] except KeyError: for key_i, val_i in value: kwargs[key_i] = val_i - if hasattr(self, 'im'): + if hasattr(self, "im"): self.im.remove() - self.im = self.ax.imshow(data*convert_flux, **kwargs) - if hasattr(self, 'cr'): + self.im = self.ax.imshow(data * convert_flux, **kwargs) + if hasattr(self, "cr"): self.cr[0].set_data(*wcs.wcs.crpix) else: - self.cr = self.ax.plot(*wcs.wcs.crpix, 'r+') + self.cr = self.ax.plot(*wcs.wcs.crpix, "r+") self.fig.canvas.draw_idle() return self.im @@ -1261,20 +1755,19 @@ class crop_map(object): def crpix_in_RS(self): crpix = self.wcs.wcs.crpix x_lim, y_lim = self.RSextent[:2], self.RSextent[2:] - if (crpix[0] > x_lim[0] and crpix[0] < x_lim[1]): - if (crpix[1] > y_lim[0] and crpix[1] < y_lim[1]): + if crpix[0] > x_lim[0] and crpix[0] < x_lim[1]: + if crpix[1] > y_lim[0] and crpix[1] < y_lim[1]: return True return False def reset_crop(self, event): self.ax.reset_wcs(self.wcs) - if hasattr(self, 'hdul_crop'): + if hasattr(self, "hdul_crop"): del self.hdul_crop, self.data_crop self.display() - if self.fig.canvas.manager.toolbar.mode == '': - self.rect_selector = RectangleSelector(self.ax, self.onselect_crop, - button=[1]) + if self.fig.canvas.manager.toolbar.mode == "": + self.rect_selector = RectangleSelector(self.ax, self.onselect_crop, button=[1]) self.RSextent = deepcopy(self.extent) self.RScenter = deepcopy(self.center) @@ -1290,7 +1783,7 @@ class crop_map(object): self.apply_crop(erelease) def apply_crop(self, event): - if hasattr(self, 'hdul_crop'): + if hasattr(self, "hdul_crop"): header = self.header_crop data = self.data_crop wcs = self.wcs_crop @@ -1306,17 +1799,17 @@ class crop_map(object): shape_im = extent[1::2] - extent[0::2] if (shape_im.astype(int) != shape).any() and (self.RSextent != self.extent).any(): # Update WCS and header in new cropped image - crpix = np.array(wcs.wcs.crpix) + # crpix = np.array(wcs.wcs.crpix) self.wcs_crop = wcs.deepcopy() self.wcs_crop.array_shape = shape if self.crpix_in_RS: self.wcs_crop.wcs.crpix = np.array(self.wcs_crop.wcs.crpix) - self.RSextent[::2] else: self.wcs_crop.wcs.crval = wcs.wcs_pix2world([self.RScenter], 1)[0] - self.wcs_crop.wcs.crpix = self.RScenter-self.RSextent[::2] + self.wcs_crop.wcs.crpix = self.RScenter - self.RSextent[::2] # Crop dataset - self.data_crop = deepcopy(data[vertex[2]:vertex[3], vertex[0]:vertex[1]]) + self.data_crop = deepcopy(data[vertex[2] : vertex[3], vertex[0] : vertex[1]]) # Write cropped map to new HDUList self.header_crop = deepcopy(header) @@ -1327,29 +1820,27 @@ class crop_map(object): self.ax.reset_wcs(self.wcs_crop) self.display(data=self.data_crop, wcs=self.wcs_crop) - xlim, ylim = self.RSextent[1::2]-self.RSextent[0::2] + xlim, ylim = self.RSextent[1::2] - self.RSextent[0::2] self.ax.set_xlim(0, xlim) self.ax.set_ylim(0, ylim) - if self.fig.canvas.manager.toolbar.mode == '': - self.rect_selector = RectangleSelector(self.ax, self.onselect_crop, - button=[1]) + if self.fig.canvas.manager.toolbar.mode == "": + self.rect_selector = RectangleSelector(self.ax, self.onselect_crop, button=[1]) self.fig.canvas.draw_idle() def on_close(self, event) -> None: - if not hasattr(self, 'hdul_crop'): + if not hasattr(self, "hdul_crop"): self.hdul_crop = self.hdul self.rect_selector.disconnect_events() self.cropped = True def crop(self) -> None: - if self.fig.canvas.manager.toolbar.mode == '': - self.rect_selector = RectangleSelector(self.ax, self.onselect_crop, - button=[1]) + if self.fig.canvas.manager.toolbar.mode == "": + self.rect_selector = RectangleSelector(self.ax, self.onselect_crop, button=[1]) self.bapply.on_clicked(self.apply_crop) self.breset.on_clicked(self.reset_crop) - self.fig.canvas.mpl_connect('close_event', self.on_close) + self.fig.canvas.mpl_connect("close_event", self.on_close) plt.show() def write_to(self, filename): @@ -1366,7 +1857,7 @@ class crop_Stokes(crop_map): """ Redefine apply_crop method for the Stokes HDUList. """ - if hasattr(self, 'hdul_crop'): + if hasattr(self, "hdul_crop"): hdul = self.hdul_crop data = self.data_crop wcs = self.wcs_crop @@ -1383,25 +1874,26 @@ class crop_Stokes(crop_map): if (shape_im.astype(int) != shape).any() and (self.RSextent != self.extent).any(): # Update WCS and header in new cropped image self.hdul_crop = deepcopy(hdul) - crpix = np.array(wcs.wcs.crpix) + self.data_crop = deepcopy(data) + # crpix = np.array(wcs.wcs.crpix) self.wcs_crop = wcs.deepcopy() self.wcs_crop.array_shape = shape if self.crpix_in_RS: self.wcs_crop.wcs.crpix = np.array(self.wcs_crop.wcs.crpix) - self.RSextent[::2] else: self.wcs_crop.wcs.crval = wcs.wcs_pix2world([self.RScenter], 1)[0] - self.wcs_crop.wcs.crpix = self.RScenter-self.RSextent[::2] + self.wcs_crop.wcs.crpix = self.RScenter - self.RSextent[::2] # Crop dataset for dataset in self.hdul_crop: - if dataset.header['datatype'] == 'IQU_cov_matrix': + if dataset.header["datatype"] == "IQU_cov_matrix": stokes_cov = np.zeros((3, 3, shape[1], shape[0])) for i in range(3): for j in range(3): - stokes_cov[i, j] = deepcopy(dataset.data[i, j][vertex[2]:vertex[3], vertex[0]:vertex[1]]) + stokes_cov[i, j] = deepcopy(dataset.data[i, j][vertex[2] : vertex[3], vertex[0] : vertex[1]]) dataset.data = stokes_cov else: - dataset.data = deepcopy(dataset.data[vertex[2]:vertex[3], vertex[0]:vertex[1]]) + dataset.data = deepcopy(dataset.data[vertex[2] : vertex[3], vertex[0] : vertex[1]]) dataset.header.update(self.wcs_crop.to_header()) self.data_crop = self.hdul_crop[0].data @@ -1410,45 +1902,49 @@ class crop_Stokes(crop_map): self.ax.reset_wcs(self.wcs_crop) self.display(data=self.data_crop, wcs=self.wcs_crop) - xlim, ylim = self.RSextent[1::2]-self.RSextent[0::2] + xlim, ylim = self.RSextent[1::2] - self.RSextent[0::2] self.ax.set_xlim(0, xlim) self.ax.set_ylim(0, ylim) else: self.on_close(event) - if self.fig.canvas.manager.toolbar.mode == '': - self.rect_selector = RectangleSelector(self.ax, self.onselect_crop, - button=[1]) + if self.fig.canvas.manager.toolbar.mode == "": + self.rect_selector = RectangleSelector(self.ax, self.onselect_crop, button=[1]) # Update integrated values - mask = np.logical_and(self.hdul_crop['data_mask'].data.astype(bool), self.hdul_crop[0].data > 0) - I_diluted = self.hdul_crop['i_stokes'].data[mask].sum() - Q_diluted = self.hdul_crop['q_stokes'].data[mask].sum() - U_diluted = self.hdul_crop['u_stokes'].data[mask].sum() - I_diluted_err = np.sqrt(np.sum(self.hdul_crop['iqu_cov_matrix'].data[0, 0][mask])) - Q_diluted_err = np.sqrt(np.sum(self.hdul_crop['iqu_cov_matrix'].data[1, 1][mask])) - U_diluted_err = np.sqrt(np.sum(self.hdul_crop['iqu_cov_matrix'].data[2, 2][mask])) - IQ_diluted_err = np.sqrt(np.sum(self.hdul_crop['iqu_cov_matrix'].data[0, 1][mask]**2)) - IU_diluted_err = np.sqrt(np.sum(self.hdul_crop['iqu_cov_matrix'].data[0, 2][mask]**2)) - QU_diluted_err = np.sqrt(np.sum(self.hdul_crop['iqu_cov_matrix'].data[1, 2][mask]**2)) + mask = np.logical_and(self.hdul_crop["data_mask"].data.astype(bool), self.hdul_crop[0].data > 0) + I_diluted = self.hdul_crop["i_stokes"].data[mask].sum() + Q_diluted = self.hdul_crop["q_stokes"].data[mask].sum() + U_diluted = self.hdul_crop["u_stokes"].data[mask].sum() + I_diluted_err = np.sqrt(np.sum(self.hdul_crop["iqu_cov_matrix"].data[0, 0][mask])) + Q_diluted_err = np.sqrt(np.sum(self.hdul_crop["iqu_cov_matrix"].data[1, 1][mask])) + U_diluted_err = np.sqrt(np.sum(self.hdul_crop["iqu_cov_matrix"].data[2, 2][mask])) + IQ_diluted_err = np.sqrt(np.sum(self.hdul_crop["iqu_cov_matrix"].data[0, 1][mask] ** 2)) + IU_diluted_err = np.sqrt(np.sum(self.hdul_crop["iqu_cov_matrix"].data[0, 2][mask] ** 2)) + QU_diluted_err = np.sqrt(np.sum(self.hdul_crop["iqu_cov_matrix"].data[1, 2][mask] ** 2)) - P_diluted = np.sqrt(Q_diluted**2+U_diluted**2)/I_diluted - P_diluted_err = (1./I_diluted)*np.sqrt((Q_diluted**2*Q_diluted_err**2 + U_diluted**2*U_diluted_err**2 + 2.*Q_diluted*U_diluted*QU_diluted_err)/(Q_diluted**2 + U_diluted ** - 2) + ((Q_diluted/I_diluted)**2 + (U_diluted/I_diluted)**2)*I_diluted_err**2 - 2.*(Q_diluted/I_diluted)*IQ_diluted_err - 2.*(U_diluted/I_diluted)*IU_diluted_err) + P_diluted = np.sqrt(Q_diluted**2 + U_diluted**2) / I_diluted + P_diluted_err = (1.0 / I_diluted) * np.sqrt( + (Q_diluted**2 * Q_diluted_err**2 + U_diluted**2 * U_diluted_err**2 + 2.0 * Q_diluted * U_diluted * QU_diluted_err) / (Q_diluted**2 + U_diluted**2) + + ((Q_diluted / I_diluted) ** 2 + (U_diluted / I_diluted) ** 2) * I_diluted_err**2 + - 2.0 * (Q_diluted / I_diluted) * IQ_diluted_err + - 2.0 * (U_diluted / I_diluted) * IU_diluted_err + ) - PA_diluted = princ_angle((90./np.pi)*np.arctan2(U_diluted, Q_diluted)) - PA_diluted_err = (90./(np.pi*(Q_diluted**2 + U_diluted**2)))*np.sqrt(U_diluted**2*Q_diluted_err ** - 2 + Q_diluted**2*U_diluted_err**2 - 2.*Q_diluted*U_diluted*QU_diluted_err) + PA_diluted = princ_angle((90.0 / np.pi) * np.arctan2(U_diluted, Q_diluted)) + PA_diluted_err = (90.0 / (np.pi * (Q_diluted**2 + U_diluted**2))) * np.sqrt( + U_diluted**2 * Q_diluted_err**2 + Q_diluted**2 * U_diluted_err**2 - 2.0 * Q_diluted * U_diluted * QU_diluted_err + ) for dataset in self.hdul_crop: - dataset.header['P_int'] = (P_diluted, 'Integrated polarization degree') - dataset.header['P_int_err'] = (np.ceil(P_diluted_err*1000.)/1000., 'Integrated polarization degree error') - dataset.header['PA_int'] = (PA_diluted, 'Integrated polarization angle') - dataset.header['PA_int_err'] = (np.ceil(PA_diluted_err*10.)/10., 'Integrated polarization angle error') + dataset.header["P_int"] = (P_diluted, "Integrated polarization degree") + dataset.header["P_int_err"] = (np.ceil(P_diluted_err * 1000.0) / 1000.0, "Integrated polarization degree error") + dataset.header["PA_int"] = (PA_diluted, "Integrated polarization angle") + dataset.header["PA_int_err"] = (np.ceil(PA_diluted_err * 10.0) / 10.0, "Integrated polarization angle error") self.fig.canvas.draw_idle() @property def data_mask(self): - return self.hdul_crop['data_mask'].data.astype(int) + return self.hdul_crop["data_mask"].data.astype(int) class image_lasso_selector(object): @@ -1458,7 +1954,7 @@ class image_lasso_selector(object): """ self.selected = False self.img = img - self.vmin, self.vmax = 0., np.max(self.img[self.img > 0.]) + self.vmin, self.vmax = 0.0, np.max(self.img[self.img > 0.0]) plt.ioff() # see https://github.com/matplotlib/matplotlib/issues/17013 if fig is None: self.fig = plt.figure(figsize=(15, 15)) @@ -1466,16 +1962,16 @@ class image_lasso_selector(object): self.fig = fig if ax is None: self.ax = self.fig.gca() - self.mask_alpha = 1. + self.mask_alpha = 1.0 self.embedded = False else: 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} + lineprops = {"color": "grey", "linewidth": 1, "alpha": 0.8} self.lasso = LassoSelector(self.ax, self.onselect, props=lineprops, useblit=False) self.lasso.set_visible(True) @@ -1484,11 +1980,11 @@ class image_lasso_selector(object): xv, yv = np.meshgrid(pix_y, pix_x) self.pix = np.vstack((xv.flatten(), yv.flatten())).T - self.fig.canvas.mpl_connect('close_event', self.on_close) + self.fig.canvas.mpl_connect("close_event", self.on_close) plt.show() def on_close(self, event=None) -> None: - if not hasattr(self, 'mask'): + if not hasattr(self, "mask"): self.mask = np.zeros(self.img.shape[:2], dtype=bool) self.lasso.disconnect_events() self.selected = True @@ -1501,15 +1997,15 @@ 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) self.mask[self.indices] = True - if hasattr(self, 'cont'): + if hasattr(self, "cont"): for coll in self.cont.collections: coll.remove() - self.cont = self.ax.contour(self.mask.astype(float), levels=[0.5], colors='white', linewidths=1) + self.cont = self.ax.contour(self.mask.astype(float), levels=[0.5], colors="white", linewidths=1) if not self.embedded: self.displayed.set_data(array) self.fig.canvas.draw_idle() @@ -1518,13 +2014,13 @@ class image_lasso_selector(object): class slit(object): - def __init__(self, img, cdelt=np.array([1., 1.]), width=1., height=2., angle=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) """ self.selected = False self.img = img - self.vmin, self.vmax = 0., np.max(self.img[self.img > 0.]) + self.vmin, self.vmax = 0.0, np.max(self.img[self.img > 0.0]) plt.ioff() # see https://github.com/matplotlib/matplotlib/issues/17013 if fig is None: self.fig = plt.figure(figsize=(15, 15)) @@ -1532,40 +2028,40 @@ class slit(object): self.fig = fig if ax is None: self.ax = self.fig.gca() - self.mask_alpha = 1. + self.mask_alpha = 1.0 self.embedded = False else: 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() xx, yy = np.indices(self.img.shape) self.pix = np.vstack((xx.flatten(), yy.flatten())).T - self.x0, self.y0 = np.array(self.img.shape)/2. + self.x0, self.y0 = np.array(self.img.shape) / 2.0 self.cdelt = cdelt - self.width = width/np.abs(self.cdelt).max()/3600. - self.height = height/np.abs(self.cdelt).max()/3600. + self.width = width / np.abs(self.cdelt).max() / 3600.0 + self.height = height / np.abs(self.cdelt).max() / 3600.0 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_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.ax.add_patch(self.rect) - self.fig.canvas.mpl_connect('button_press_event', self.on_press) - self.fig.canvas.mpl_connect('button_release_event', self.on_release) - self.fig.canvas.mpl_connect('motion_notify_event', self.on_move) - self.fig.canvas.mpl_connect('close_event', self.on_close) + self.fig.canvas.mpl_connect("button_press_event", self.on_press) + self.fig.canvas.mpl_connect("button_release_event", self.on_release) + self.fig.canvas.mpl_connect("motion_notify_event", self.on_move) + self.fig.canvas.mpl_connect("close_event", self.on_close) self.x0, self.y0 = self.rect.xy self.pressevent = None plt.show() def on_close(self, event=None) -> None: - if not hasattr(self, 'mask'): + if not hasattr(self, "mask"): self.mask = np.zeros(self.img.shape[:2], dtype=bool) self.selected = True @@ -1593,12 +2089,12 @@ class slit(object): self.fig.canvas.draw_idle() def update_width(self, width): - self.width = width/np.abs(self.cdelt).max()/3600 + self.width = width / np.abs(self.cdelt).max() / 3600 self.rect.set_width(self.width) self.fig.canvas.draw_idle() def update_height(self, height): - self.height = height/np.abs(self.cdelt).max()/3600 + self.height = height / np.abs(self.cdelt).max() / 3600 self.rect.set_height(self.height) self.fig.canvas.draw_idle() @@ -1608,24 +2104,24 @@ class slit(object): self.fig.canvas.draw_idle() def update_mask(self): - if hasattr(self, 'displayed'): + if hasattr(self, "displayed"): try: 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) for p in self.pix: - self.mask[tuple(p)] = (np.abs(np.dot(rot2D(-self.angle), p-self.rect.get_center()[::-1])) < (self.height/2., self.width/2.)).all() - if hasattr(self, 'cont'): + self.mask[tuple(p)] = (np.abs(np.dot(rot2D(-self.angle), p - self.rect.get_center()[::-1])) < (self.height / 2.0, self.width / 2.0)).all() + if hasattr(self, "cont"): for coll in self.cont.collections: try: coll.remove() except AttributeError: return - self.cont = self.ax.contour(self.mask.astype(float), levels=[0.5], colors='white', linewidths=1) + self.cont = self.ax.contour(self.mask.astype(float), levels=[0.5], colors="white", linewidths=1) if not self.embedded: self.displayed.set_data(array) self.fig.canvas.draw_idle() @@ -1634,13 +2130,13 @@ class slit(object): class aperture(object): - def __init__(self, img, cdelt=np.array([1., 1.]), radius=1., fig=None, ax=None): + def __init__(self, img, cdelt=np.array([1.0, 1.0]), radius=1.0, fig=None, ax=None): """ img must have shape (X, Y) """ self.selected = False self.img = img - self.vmin, self.vmax = 0., np.max(self.img[self.img > 0.]) + self.vmin, self.vmax = 0.0, np.max(self.img[self.img > 0.0]) plt.ioff() # see https://github.com/matplotlib/matplotlib/issues/17013 if fig is None: self.fig = plt.figure(figsize=(15, 15)) @@ -1648,37 +2144,37 @@ class aperture(object): self.fig = fig if ax is None: self.ax = self.fig.gca() - self.mask_alpha = 1. + self.mask_alpha = 1.0 self.embedded = False else: 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() xx, yy = np.indices(self.img.shape) self.pix = np.vstack((xx.flatten(), yy.flatten())).T - self.x0, self.y0 = np.array(self.img.shape)/2. - if np.abs(cdelt).max() != 1.: + self.x0, self.y0 = np.array(self.img.shape) / 2.0 + if np.abs(cdelt).max() != 1.0: self.cdelt = cdelt - self.radius = radius/np.abs(self.cdelt).max()/3600. + self.radius = radius / np.abs(self.cdelt).max() / 3600.0 - self.circ = Circle((self.x0, self.y0), self.radius, alpha=0.8, ec='grey', fc='none') + self.circ = Circle((self.x0, self.y0), self.radius, alpha=0.8, ec="grey", fc="none") self.ax.add_patch(self.circ) - self.fig.canvas.mpl_connect('button_press_event', self.on_press) - self.fig.canvas.mpl_connect('button_release_event', self.on_release) - self.fig.canvas.mpl_connect('motion_notify_event', self.on_move) - self.fig.canvas.mpl_connect('close_event', self.on_close) + self.fig.canvas.mpl_connect("button_press_event", self.on_press) + self.fig.canvas.mpl_connect("button_release_event", self.on_release) + self.fig.canvas.mpl_connect("motion_notify_event", self.on_move) + self.fig.canvas.mpl_connect("close_event", self.on_close) self.x0, self.y0 = self.circ.center self.pressevent = None plt.show() def on_close(self, event=None) -> None: - if not hasattr(self, 'mask'): + if not hasattr(self, "mask"): self.mask = np.zeros(self.img.shape[:2], dtype=bool) self.selected = True @@ -1706,29 +2202,29 @@ class aperture(object): self.fig.canvas.draw_idle() def update_radius(self, radius): - self.radius = radius/np.abs(self.cdelt).max()/3600 + self.radius = radius / np.abs(self.cdelt).max() / 3600 self.circ.set_radius(self.radius) self.fig.canvas.draw_idle() def update_mask(self): - if hasattr(self, 'displayed'): + if hasattr(self, "displayed"): try: 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]) x0, y0 = self.circ.center - self.mask = np.sqrt((xx-x0)**2+(yy-y0)**2) < self.radius - if hasattr(self, 'cont'): + self.mask = np.sqrt((xx - x0) ** 2 + (yy - y0) ** 2) < self.radius + if hasattr(self, "cont"): for coll in self.cont.collections: try: coll.remove() except AttributeError: return - self.cont = self.ax.contour(self.mask.astype(float), levels=[0.5], colors='white', linewidths=1) + self.cont = self.ax.contour(self.mask.astype(float), levels=[0.5], colors="white", linewidths=1) if not self.embedded: self.displayed.set_data(array) self.fig.canvas.draw_idle() @@ -1741,8 +2237,7 @@ class pol_map(object): Class to interactively study polarization maps. """ - def __init__(self, Stokes, SNRp_cut=3., SNRi_cut=3., flux_lim=None, selection=None): - + def __init__(self, Stokes, SNRp_cut=3.0, SNRi_cut=3.0, step_vec=1, scale_vec=3.0, flux_lim=None, selection=None): if isinstance(Stokes, str): Stokes = fits.open(Stokes) self.Stokes = deepcopy(Stokes) @@ -1754,15 +2249,16 @@ class pol_map(object): self.region = None self.data = None self.display_selection = selection - self.vec_scale = 2. + self.step_vec = step_vec + self.scale_vec = scale_vec # 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[0].header["targname"] + self.pivot_wav = self.Stokes[0].header["photplam"] + self.map_convert = self.Stokes[0].header["photflam"] # Create figure - plt.rcParams.update({'font.size': 10}) + plt.rcParams.update({"font.size": 10}) self.fig, self.ax = plt.subplots(figsize=(10, 10), subplot_kw=dict(projection=self.wcs)) self.fig.subplots_adjust(hspace=0, wspace=0, right=1.02) self.ax_cosmetics() @@ -1779,11 +2275,11 @@ class pol_map(object): ax_P_cut = self.fig.add_axes([0.120, 0.055, 0.230, 0.01]) ax_vec_sc = self.fig.add_axes([0.240, 0.030, 0.110, 0.01]) ax_snr_reset = self.fig.add_axes([0.080, 0.020, 0.05, 0.02]) - SNRi_max = np.max(self.I[self.IQU_cov[0, 0] > 0.]/np.sqrt(self.IQU_cov[0, 0][self.IQU_cov[0, 0] > 0.])) - SNRp_max = np.max(self.P[self.s_P > 0.]/self.s_P[self.s_P > 0.]) - s_I_cut = Slider(ax_I_cut, r"$SNR^{I}_{cut}$", 1., int(SNRi_max*0.95), valstep=1, valinit=self.SNRi_cut) - s_P_cut = Slider(ax_P_cut, r"$SNR^{P}_{cut}$", 1., int(SNRp_max*0.95), valstep=1, valinit=self.SNRp_cut) - s_vec_sc = Slider(ax_vec_sc, r"Vectors scale", 1., 6., valstep=1, valinit=self.vec_scale) + 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) + s_P_cut = Slider(ax_P_cut, r"$SNR^{P}_{cut}$", 1.0, int(SNRp_max * 0.95), valstep=1, valinit=self.SNRp_cut) + s_vec_sc = Slider(ax_vec_sc, r"Vectors 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) @@ -1800,7 +2296,7 @@ class pol_map(object): self.fig.canvas.draw_idle() def update_vecsc(val): - self.vec_scale = val + self.scale_vec = val self.pol_vector() self.ax_cosmetics() self.fig.canvas.draw_idle() @@ -1865,7 +2361,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.) + 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: @@ -1889,8 +2385,8 @@ class pol_map(object): self.fig.canvas.draw_idle() def update_aperture(val): - if hasattr(self, 'select_instance'): - if hasattr(self.select_instance, 'radius'): + if hasattr(self, "select_instance"): + if hasattr(self.select_instance, "radius"): self.select_instance.update_radius(val) else: self.selected = True @@ -1921,9 +2417,9 @@ 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., valstep=1e-2, valinit=1.) - s_slit_height = Slider(ax_slit_height, r"$H_{slit}$", np.ceil(self.wcs.wcs.cdelt.max()/1.33*3.6e5)/1e2, 7., valstep=1e-2, valinit=1.) - s_slit_angle = Slider(ax_slit_angle, r"$\theta_{slit}$", 0., 90., valstep=1., valinit=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): if self.data is None: @@ -1941,52 +2437,59 @@ class pol_map(object): else: 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.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.select_instance.rect.set_visible(True) self.fig.canvas.draw_idle() def update_slit_w(val): - if hasattr(self, 'select_instance'): - if hasattr(self.select_instance, 'width'): + if hasattr(self, "select_instance"): + if hasattr(self.select_instance, "width"): self.select_instance.update_width(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.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 + ) 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.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.fig.canvas.draw_idle() def update_slit_h(val): - if hasattr(self, 'select_instance'): - if hasattr(self.select_instance, 'height'): + if hasattr(self, "select_instance"): + if hasattr(self.select_instance, "height"): self.select_instance.update_height(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.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 + ) 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.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.fig.canvas.draw_idle() def update_slit_a(val): - if hasattr(self, 'select_instance'): - if hasattr(self.select_instance, 'angle'): + if hasattr(self, "select_instance"): + if hasattr(self.select_instance, "angle"): self.select_instance.update_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.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 + ) 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.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.fig.canvas.draw_idle() def reset_slit(event): @@ -2055,7 +2558,7 @@ class pol_map(object): b_save = Button(ax_save, "Save") b_save.label.set_fontsize(8) ax_text_save = self.fig.add_axes([0.3, 0.020, 0.5, 0.025], visible=False) - text_save = TextBox(ax_text_save, "Save to:", initial='') + text_save = TextBox(ax_text_save, "Save to:", initial="") def saveplot(event): ax_text_save.set(visible=True) @@ -2069,18 +2572,18 @@ 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)) + if expression != "": + 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) self.pol_int(fig=save_fig, ax=save_ax) save_fig.suptitle(r"{0:s} with $SNR_{{p}} \geq$ {1:d} and $SNR_{{I}} \geq$ {2:d}".format(self.targ, int(self.SNRp), int(self.SNRi))) - if not expression[-4:] in ['.png', '.jpg', '.pdf']: - expression += '.pdf' - save_fig.savefig(expression, bbox_inches='tight', dpi=200) + if not expression[-4:] in [".png", ".jpg", ".pdf"]: + expression += ".pdf" + save_fig.savefig(expression, bbox_inches="tight", dpi=200) plt.close(save_fig) - text_save.set_val('') + text_save.set_val("") ax_snr_reset.set(visible=True) ax_vec_sc.set(visible=True) ax_save.set(visible=True) @@ -2094,7 +2597,7 @@ class pol_map(object): b_dump = Button(ax_dump, "Dump") b_dump.label.set_fontsize(8) ax_text_dump = self.fig.add_axes([0.3, 0.020, 0.5, 0.025], visible=False) - text_dump = TextBox(ax_text_dump, "Dump to:", initial='') + text_dump = TextBox(ax_text_dump, "Dump to:", initial="") def dump(event): ax_text_dump.set(visible=True) @@ -2105,10 +2608,10 @@ class pol_map(object): self.fig.canvas.draw_idle() shape = np.array(self.I.shape) - center = (shape/2).astype(int) - cdelt_arcsec = self.wcs.wcs.cdelt*3600 + 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] @@ -2116,19 +2619,20 @@ class pol_map(object): dump_list = [] 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]]) + 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]] + ) self.data_dump = np.array(dump_list) b_dump.on_clicked(dump) def submit_dump(expression): ax_text_dump.set(visible=False) - if expression != '': - if not expression[-4:] in ['.txt', '.dat']: - expression += '.txt' + if expression != "": + if not expression[-4:] in [".txt", ".dat"]: + expression += ".txt" np.savetxt(expression, self.data_dump) - text_dump.set_val('') + text_dump.set_val("") ax_snr_reset.set(visible=True) ax_vec_sc.set(visible=True) ax_save.set(visible=True) @@ -2152,39 +2656,45 @@ class pol_map(object): b_snrp = Button(ax_snrp, r"$P / \sigma_{P}$") def d_tf(event): - self.display_selection = 'total_flux' + self.display_selection = "total_flux" self.display() self.pol_int() + b_tf.on_clicked(d_tf) def d_pf(event): - self.display_selection = 'pol_flux' + self.display_selection = "pol_flux" self.display() self.pol_int() + b_pf.on_clicked(d_pf) def d_p(event): - self.display_selection = 'pol_deg' + self.display_selection = "pol_deg" self.display() self.pol_int() + b_p.on_clicked(d_p) def d_pa(event): - self.display_selection = 'pol_ang' + self.display_selection = "pol_ang" self.display() self.pol_int() + b_pa.on_clicked(d_pa) def d_snri(event): - self.display_selection = 'snri' + self.display_selection = "snri" self.display() self.pol_int() + b_snri.on_clicked(d_snri) def d_snrp(event): - self.display_selection = 'snrp' + self.display_selection = "snrp" self.display() self.pol_int() + b_snrp.on_clicked(d_snrp) plt.show() @@ -2195,75 +2705,107 @@ class pol_map(object): @property def I(self): - return self.Stokes['I_STOKES'].data + return self.Stokes["I_STOKES"].data @property def Q(self): - return self.Stokes['Q_STOKES'].data + return self.Stokes["Q_STOKES"].data @property def U(self): - return self.Stokes['U_STOKES'].data + return self.Stokes["U_STOKES"].data @property def IQU_cov(self): - return self.Stokes['IQU_COV_MATRIX'].data + return self.Stokes["IQU_COV_MATRIX"].data @property def P(self): - return self.Stokes['POL_DEG_DEBIASED'].data + return self.Stokes["POL_DEG_DEBIASED"].data @property def s_P(self): - return self.Stokes['POL_DEG_ERR'].data + return self.Stokes["POL_DEG_ERR"].data @property def PA(self): - return self.Stokes['POL_ANG'].data + return self.Stokes["POL_ANG"].data @property def data_mask(self): - return self.Stokes['DATA_MASK'].data + return self.Stokes["DATA_MASK"].data 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) + self.Stokes[np.argmax([self.Stokes[i].header["datatype"] == "Data_mask" for i in range(len(self.Stokes))])].data = mask.astype(float) @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[self.s_P > 0.] = self.P[self.s_P > 0.] / self.s_P[self.s_P > 0.] > self.SNRp - SNRi_mask[s_I > 0.] = self.I[s_I > 0.] / s_I[s_I > 0.] > self.SNRi + SNRp_mask[self.s_P > 0.0] = self.P[self.s_P > 0.0] / self.s_P[self.s_P > 0.0] > self.SNRp + SNRi_mask[s_I > 0.0] = self.I[s_I > 0.0] / s_I[s_I > 0.0] > self.SNRi return np.logical_and(SNRi_mask, SNRp_mask) def ax_cosmetics(self, ax=None): if ax is None: ax = self.ax - ax.set(aspect='equal', fc='black') + ax.set(aspect="equal", fc="black") - ax.coords.grid(True, color='white', ls='dotted', alpha=0.5) - ax.coords[0].set_axislabel('Right Ascension (J2000)') - ax.coords[0].set_axislabel_position('t') - ax.coords[0].set_ticklabel_position('t') - ax.set_ylabel('Declination (J2000)', labelpad=-1) + ax.coords.grid(True, color="white", ls="dotted", alpha=0.5) + ax.coords[0].set_axislabel("Right Ascension (J2000)") + ax.coords[0].set_axislabel_position("t") + ax.coords[0].set_ticklabel_position("t") + ax.set_ylabel("Declination (J2000)", labelpad=-1) # Display scales and orientation fontprops = fm.FontProperties(size=14) - px_size = self.wcs.wcs.cdelt[0]*3600. - if hasattr(self, 'px_sc'): + px_size = self.wcs.wcs.cdelt[0] * 3600.0 + if hasattr(self, "px_sc"): self.px_sc.remove() - self.px_sc = AnchoredSizeBar(ax.transData, 1./px_size, '1 arcsec', 3, pad=0.5, sep=5, borderpad=0.5, - frameon=False, size_vertical=0.005, color='white', fontproperties=fontprops) + self.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="white", + fontproperties=fontprops, + ) ax.add_artist(self.px_sc) - if hasattr(self, 'pol_sc'): + if hasattr(self, "pol_sc"): self.pol_sc.remove() - self.pol_sc = AnchoredSizeBar(ax.transData, self.vec_scale, r"$P$= 100%", 4, pad=0.5, sep=5, borderpad=0.5, - frameon=False, size_vertical=0.005, color='white', fontproperties=fontprops) + if self.scale_vec != 0: + scale_vec = self.scale_vec + 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.add_artist(self.pol_sc) - if hasattr(self, 'north_dir'): + if hasattr(self, "north_dir"): self.north_dir.remove() - self.north_dir = AnchoredDirectionArrows(ax.transAxes, "E", "N", length=-0.08, fontsize=0.025, loc=1, aspect_ratio=-1, sep_y=0.01, sep_x=0.01, back_length=0., head_length=10., head_width=10., - 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}) + self.north_dir = AnchoredDirectionArrows( + ax.transAxes, + "E", + "N", + length=-0.08, + fontsize=0.025, + loc=1, + aspect_ratio=-1, + sep_y=0.01, + sep_x=0.01, + back_length=0.0, + head_length=10.0, + head_width=10.0, + 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}, + ) ax.add_artist(self.north_dir) def display(self, fig=None, ax=None, flux_lim=None): @@ -2272,93 +2814,128 @@ class pol_map(object): self.display_selection = "total_flux" 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.lower() in ["total_flux"]: + self.data = self.I * self.map_convert if flux_lim is None: - vmin, vmax = 1./2.*np.median(self.data[self.data > 0.]), np.max(self.data[self.data > 0.]) + 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) 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 + 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./2.*np.median(self.I[self.I > 0.]*self.map_convert), np.max(self.I[self.I > 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 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. - vmin, vmax = 0., np.max(self.data[self.P > self.s_P]) + 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]) label = r"$P$ [%]" - elif self.display_selection.lower() in ['pol_ang']: + elif self.display_selection.lower() in ["pol_ang"]: self.data = princ_angle(self.PA) - vmin, vmax = 0, 180. + vmin, vmax = 0, 180.0 label = r"$\theta_{P}$ [°]" - elif self.display_selection.lower() in ['snri']: + 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.] = self.I[s_I > 0.]/s_I[s_I > 0.] + SNRi[s_I > 0.0] = self.I[s_I > 0.0] / s_I[s_I > 0.0] self.data = SNRi - vmin, vmax = 0., np.max(self.data[self.data > 0.]) + vmin, vmax = 0.0, np.max(self.data[self.data > 0.0]) label = r"$I_{Stokes}/\sigma_{I}$" - elif self.display_selection.lower() in ['snrp']: + elif self.display_selection.lower() in ["snrp"]: SNRp = np.zeros(self.P.shape) - SNRp[self.s_P > 0.] = self.P[self.s_P > 0.]/self.s_P[self.s_P > 0.] + 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., np.max(self.data[self.data > 0.]) + vmin, vmax = 0.0, np.max(self.data[self.data > 0.0]) label = r"$P/\sigma_{P}$" if fig is None: fig = self.fig if ax is None: ax = self.ax - if hasattr(self, 'cbar'): + if hasattr(self, "cbar"): self.cbar.remove() - if hasattr(self, 'im'): + if hasattr(self, "im"): self.im.remove() if norm is not None: - self.im = ax.imshow(self.data, norm=norm, aspect='equal', cmap='inferno') + 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') - plt.rcParams.update({'font.size': 14}) + self.im = ax.imshow(self.data, vmin=vmin, vmax=vmax, aspect="equal", cmap="inferno") + 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}) + 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') + 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') + 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]) - plt.rcParams.update({'font.size': 14}) + 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}) + plt.rcParams.update({"font.size": 10}) fig.canvas.draw_idle() return im def pol_vector(self, fig=None, ax=None): - P_cut = np.ones(self.P.shape)*np.nan - P_cut[self.cut] = self.P[self.cut] + P_cut = np.ones(self.P.shape) * np.nan + if self.scale_vec != 0.0: + scale_vec = self.scale_vec + P_cut[self.cut] = self.P[self.cut] + else: + 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. + self.PA*np.pi/180.), P_cut*np.sin(np.pi/2. + self.PA*np.pi/180.) + 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 if ax is None: ax = self.ax - if hasattr(self, 'quiver'): + if hasattr(self, "quiver"): self.quiver.remove() - self.quiver = ax.quiver(X, Y, XY_U, XY_V, units='xy', scale=1./self.vec_scale, scale_units='xy', pivot='mid', headwidth=0., - headlength=0., headaxislength=0., width=0.5, linewidth=0.75, color='white', edgecolor='black') + 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.4, + linewidth=0.8, + color="white", + edgecolor="black", + ) fig.canvas.draw_idle() return self.quiver else: - ax.quiver(X, Y, XY_U, XY_V, units='xy', scale=1./self.vec_scale, scale_units='xy', pivot='mid', headwidth=0., - headlength=0., headaxislength=0., width=0.5, linewidth=0.75, color='white', edgecolor='black') + 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.4, + linewidth=0.8, + color="white", + edgecolor="black", + ) fig.canvas.draw_idle() def pol_int(self, fig=None, ax=None): @@ -2366,10 +2943,10 @@ class pol_map(object): 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['P_int_err'] - PA_reg = self.Stokes[0].header['PA_int'] - PA_reg_err = self.Stokes[0].header['PA_int_err'] + P_reg = self.Stokes[0].header["P_int"] + P_reg_err = self.Stokes[0].header["P_int_err"] + PA_reg = self.Stokes[0].header["PA_int"] + PA_reg_err = self.Stokes[0].header["PA_int_err"] s_I = np.sqrt(self.IQU_cov[0, 0]) s_Q = np.sqrt(self.IQU_cov[1, 1]) @@ -2381,19 +2958,29 @@ class pol_map(object): I_cut = self.I[self.cut].sum() Q_cut = self.Q[self.cut].sum() U_cut = self.U[self.cut].sum() - I_cut_err = np.sqrt(np.sum(s_I[self.cut]**2)) - Q_cut_err = np.sqrt(np.sum(s_Q[self.cut]**2)) - U_cut_err = np.sqrt(np.sum(s_U[self.cut]**2)) - IQ_cut_err = np.sqrt(np.sum(s_IQ[self.cut]**2)) - IU_cut_err = np.sqrt(np.sum(s_IU[self.cut]**2)) - QU_cut_err = np.sqrt(np.sum(s_QU[self.cut]**2)) + I_cut_err = np.sqrt(np.sum(s_I[self.cut] ** 2)) + Q_cut_err = np.sqrt(np.sum(s_Q[self.cut] ** 2)) + U_cut_err = np.sqrt(np.sum(s_U[self.cut] ** 2)) + IQ_cut_err = np.sqrt(np.sum(s_IQ[self.cut] ** 2)) + IU_cut_err = np.sqrt(np.sum(s_IU[self.cut] ** 2)) + QU_cut_err = np.sqrt(np.sum(s_QU[self.cut] ** 2)) - P_cut = np.sqrt(Q_cut**2+U_cut**2)/I_cut - P_cut_err = np.sqrt((Q_cut**2*Q_cut_err**2 + U_cut**2*U_cut_err**2 + 2.*Q_cut*U_cut*QU_cut_err)/(Q_cut**2 + U_cut**2) + - ((Q_cut/I_cut)**2 + (U_cut/I_cut)**2)*I_cut_err**2 - 2.*(Q_cut/I_cut)*IQ_cut_err - 2.*(U_cut/I_cut)*IU_cut_err)/I_cut + with np.errstate(divide="ignore", invalid="ignore"): + P_cut = np.sqrt(Q_cut**2 + U_cut**2) / I_cut + P_cut_err = ( + np.sqrt( + (Q_cut**2 * Q_cut_err**2 + U_cut**2 * U_cut_err**2 + 2.0 * Q_cut * U_cut * QU_cut_err) / (Q_cut**2 + U_cut**2) + + ((Q_cut / I_cut) ** 2 + (U_cut / I_cut) ** 2) * I_cut_err**2 + - 2.0 * (Q_cut / I_cut) * IQ_cut_err + - 2.0 * (U_cut / I_cut) * IU_cut_err + ) + / I_cut + ) - PA_cut = princ_angle((90./np.pi)*np.arctan2(U_cut, Q_cut)) - PA_cut_err = (90./(np.pi*(Q_cut**2+U_cut**2)))*np.sqrt(U_cut**2*Q_cut_err**2 + Q_cut**2*U_cut_err**2 - 2.*Q_cut*U_cut*QU_cut_err) + PA_cut = princ_angle((90.0 / np.pi) * np.arctan2(U_cut, Q_cut)) + PA_cut_err = (90.0 / (np.pi * (Q_cut**2 + U_cut**2))) * np.sqrt( + U_cut**2 * Q_cut_err**2 + Q_cut**2 * U_cut_err**2 - 2.0 * Q_cut * U_cut * QU_cut_err + ) else: s_I = np.sqrt(self.IQU_cov[0, 0]) @@ -2406,39 +2993,59 @@ class pol_map(object): I_reg = self.I[self.region].sum() Q_reg = self.Q[self.region].sum() U_reg = self.U[self.region].sum() - I_reg_err = np.sqrt(np.sum(s_I[self.region]**2)) - Q_reg_err = np.sqrt(np.sum(s_Q[self.region]**2)) - U_reg_err = np.sqrt(np.sum(s_U[self.region]**2)) - IQ_reg_err = np.sqrt(np.sum(s_IQ[self.region]**2)) - IU_reg_err = np.sqrt(np.sum(s_IU[self.region]**2)) - QU_reg_err = np.sqrt(np.sum(s_QU[self.region]**2)) + I_reg_err = np.sqrt(np.sum(s_I[self.region] ** 2)) + Q_reg_err = np.sqrt(np.sum(s_Q[self.region] ** 2)) + U_reg_err = np.sqrt(np.sum(s_U[self.region] ** 2)) + IQ_reg_err = np.sqrt(np.sum(s_IQ[self.region] ** 2)) + IU_reg_err = np.sqrt(np.sum(s_IU[self.region] ** 2)) + QU_reg_err = np.sqrt(np.sum(s_QU[self.region] ** 2)) - P_reg = np.sqrt(Q_reg**2+U_reg**2)/I_reg - P_reg_err = np.sqrt((Q_reg**2*Q_reg_err**2 + U_reg**2*U_reg_err**2 + 2.*Q_reg*U_reg*QU_reg_err)/(Q_reg**2 + U_reg**2) + - ((Q_reg/I_reg)**2 + (U_reg/I_reg)**2)*I_reg_err**2 - 2.*(Q_reg/I_reg)*IQ_reg_err - 2.*(U_reg/I_reg)*IU_reg_err)/I_reg + with np.errstate(divide="ignore", invalid="ignore"): + P_reg = np.sqrt(Q_reg**2 + U_reg**2) / I_reg + P_reg_err = ( + np.sqrt( + (Q_reg**2 * Q_reg_err**2 + U_reg**2 * U_reg_err**2 + 2.0 * Q_reg * U_reg * QU_reg_err) / (Q_reg**2 + U_reg**2) + + ((Q_reg / I_reg) ** 2 + (U_reg / I_reg) ** 2) * I_reg_err**2 + - 2.0 * (Q_reg / I_reg) * IQ_reg_err + - 2.0 * (U_reg / I_reg) * IU_reg_err + ) + / I_reg + ) - PA_reg = princ_angle((90./np.pi)*np.arctan2(U_reg, Q_reg)) - PA_reg_err = (90./(np.pi*(Q_reg**2+U_reg**2)))*np.sqrt(U_reg**2*Q_reg_err**2 + Q_reg**2*U_reg_err**2 - 2.*Q_reg*U_reg*QU_reg_err) + PA_reg = princ_angle((90.0 / np.pi) * np.arctan2(U_reg, Q_reg)) + PA_reg_err = (90.0 / (np.pi * (Q_reg**2 + U_reg**2))) * np.sqrt( + U_reg**2 * Q_reg_err**2 + Q_reg**2 * U_reg_err**2 - 2.0 * Q_reg * U_reg * QU_reg_err + ) new_cut = np.logical_and(self.region, self.cut) I_cut = self.I[new_cut].sum() Q_cut = self.Q[new_cut].sum() U_cut = self.U[new_cut].sum() - I_cut_err = np.sqrt(np.sum(s_I[new_cut]**2)) - Q_cut_err = np.sqrt(np.sum(s_Q[new_cut]**2)) - U_cut_err = np.sqrt(np.sum(s_U[new_cut]**2)) - IQ_cut_err = np.sqrt(np.sum(s_IQ[new_cut]**2)) - IU_cut_err = np.sqrt(np.sum(s_IU[new_cut]**2)) - QU_cut_err = np.sqrt(np.sum(s_QU[new_cut]**2)) + I_cut_err = np.sqrt(np.sum(s_I[new_cut] ** 2)) + Q_cut_err = np.sqrt(np.sum(s_Q[new_cut] ** 2)) + U_cut_err = np.sqrt(np.sum(s_U[new_cut] ** 2)) + IQ_cut_err = np.sqrt(np.sum(s_IQ[new_cut] ** 2)) + IU_cut_err = np.sqrt(np.sum(s_IU[new_cut] ** 2)) + QU_cut_err = np.sqrt(np.sum(s_QU[new_cut] ** 2)) - P_cut = np.sqrt(Q_cut**2+U_cut**2)/I_cut - P_cut_err = np.sqrt((Q_cut**2*Q_cut_err**2 + U_cut**2*U_cut_err**2 + 2.*Q_cut*U_cut*QU_cut_err)/(Q_cut**2 + U_cut**2) + - ((Q_cut/I_cut)**2 + (U_cut/I_cut)**2)*I_cut_err**2 - 2.*(Q_cut/I_cut)*IQ_cut_err - 2.*(U_cut/I_cut)*IU_cut_err)/I_cut + with np.errstate(divide="ignore", invalid="ignore"): + P_cut = np.sqrt(Q_cut**2 + U_cut**2) / I_cut + P_cut_err = ( + np.sqrt( + (Q_cut**2 * Q_cut_err**2 + U_cut**2 * U_cut_err**2 + 2.0 * Q_cut * U_cut * QU_cut_err) / (Q_cut**2 + U_cut**2) + + ((Q_cut / I_cut) ** 2 + (U_cut / I_cut) ** 2) * I_cut_err**2 + - 2.0 * (Q_cut / I_cut) * IQ_cut_err + - 2.0 * (U_cut / I_cut) * IU_cut_err + ) + / I_cut + ) - PA_cut = princ_angle((90./np.pi)*np.arctan2(U_cut, Q_cut)) - PA_cut_err = (90./(np.pi*(Q_cut**2+U_cut**2)))*np.sqrt(U_cut**2*Q_cut_err**2 + Q_cut**2*U_cut_err**2 - 2.*Q_cut*U_cut*QU_cut_err) + PA_cut = princ_angle((90.0 / np.pi) * np.arctan2(U_cut, Q_cut)) + PA_cut_err = (90.0 / (np.pi * (Q_cut**2 + U_cut**2))) * np.sqrt( + U_cut**2 * Q_cut_err**2 + Q_cut**2 * U_cut_err**2 - 2.0 * Q_cut * U_cut * QU_cut_err + ) - if hasattr(self, 'cont'): + if hasattr(self, "cont"): for coll in self.cont.collections: try: coll.remove() @@ -2449,39 +3056,81 @@ class pol_map(object): fig = self.fig if ax is None: ax = self.ax - if hasattr(self, 'an_int'): + if hasattr(self, "an_int"): 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))+"\n"+r"$P^{{int}}$ = {0:.1f} $\pm$ {1:.1f} %".format(P_reg*100., np.ceil(P_reg_err*1000.)/10.)+"\n"+r"$\theta_{{P}}^{{int}}$ = {0:.1f} $\pm$ {1:.1f} °".format(PA_reg, np.ceil(PA_reg_err*10.)/10.) + 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) + ) + + "\n" + + r"$P^{{int}}$ = {0:.1f} $\pm$ {1:.1f} %".format(P_reg * 100.0, np.ceil(P_reg_err * 1000.0) / 10.0) + + "\n" + + r"$\theta_{{P}}^{{int}}$ = {0:.1f} $\pm$ {1:.1f} °".format(PA_reg, np.ceil(PA_reg_err * 10.0) / 10.0) + ) self.str_cut = "" # self.str_cut = "\n"+r"$F_{{\lambda}}^{{cut}}$({0:.0f} $\AA$) = {1} $ergs \cdot cm^{{-2}} \cdot s^{{-1}} \cdot \AA^{{-1}}$".format(self.pivot_wav, sci_not(I_cut*self.map_convert, I_cut_err*self.map_convert, 2))+"\n"+r"$P^{{cut}}$ = {0:.1f} $\pm$ {1:.1f} %".format(P_cut*100., np.ceil(P_cut_err*1000.)/10.)+"\n"+r"$\theta_{{P}}^{{cut}}$ = {0:.1f} $\pm$ {1:.1f} °".format(PA_cut, np.ceil(PA_cut_err*10.)/10.) - self.an_int = ax.annotate(self.str_int+self.str_cut, color='white', fontsize=12, xy=(0.01, 1.00), xycoords='axes fraction', path_effects=[pe.withStroke(linewidth=0.5, foreground='k')], verticalalignment='top', horizontalalignment='left') + self.an_int = ax.annotate( + self.str_int + self.str_cut, + color="white", + fontsize=12, + xy=(0.01, 1.00), + xycoords="axes fraction", + path_effects=[pe.withStroke(linewidth=0.5, foreground="k")], + verticalalignment="top", + 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))+"\n"+r"$P^{{int}}$ = {0:.1f} $\pm$ {1:.1f} %".format(P_reg*100., np.ceil(P_reg_err*1000.)/10.)+"\n"+r"$\theta_{{P}}^{{int}}$ = {0:.1f} $\pm$ {1:.1f} °".format(PA_reg, np.ceil(PA_reg_err*10.)/10.) + 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) + ) + + "\n" + + r"$P^{{int}}$ = {0:.1f} $\pm$ {1:.1f} %".format(P_reg * 100.0, np.ceil(P_reg_err * 1000.0) / 10.0) + + "\n" + + r"$\theta_{{P}}^{{int}}$ = {0:.1f} $\pm$ {1:.1f} °".format(PA_reg, np.ceil(PA_reg_err * 10.0) / 10.0) + ) str_cut = "" # str_cut = "\n"+r"$F_{{\lambda}}^{{cut}}$({0:.0f} $\AA$) = {1} $ergs \cdot cm^{{-2}} \cdot s^{{-1}} \cdot \AA^{{-1}}$".format(self.pivot_wav, sci_not(I_cut*self.map_convert, I_cut_err*self.map_convert, 2))+"\n"+r"$P^{{cut}}$ = {0:.1f} $\pm$ {1:.1f} %".format(P_cut*100., np.ceil(P_cut_err*1000.)/10.)+"\n"+r"$\theta_{{P}}^{{cut}}$ = {0:.1f} $\pm$ {1:.1f} °".format(PA_cut, np.ceil(PA_cut_err*10.)/10.) - ax.annotate(str_int+str_cut, color='white', fontsize=12, xy=(0.01, 1.00), xycoords='axes fraction', path_effects=[pe.withStroke(linewidth=0.5, foreground='k')], verticalalignment='top', horizontalalignment='left') + ax.annotate( + str_int + str_cut, + color="white", + fontsize=12, + xy=(0.01, 1.00), + xycoords="axes fraction", + path_effects=[pe.withStroke(linewidth=0.5, foreground="k")], + verticalalignment="top", + 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() 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', '--snrp', metavar='snrp_cut', required=False, help='the cut in signal-to-noise for the polarization degree', type=float, default=3.) - parser.add_argument('-i', '--snri', metavar='snri_cut', required=False, help='the cut in signal-to-noise for the intensity', type=float, default=3.) - parser.add_argument('-l', '--lim', metavar='flux_lim', nargs=2, required=False, help='limits for the intensity map', default=None) + 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", "--snrp", metavar="snrp_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( + "-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("-l", "--lim", metavar="flux_lim", nargs=2, required=False, help="limits for the intensity map", default=None) args = parser.parse_args() if args.file is not None: - Stokes_UV = fits.open(args.file, mode='readonly') - p = pol_map(Stokes_UV, SNRp_cut=args.snrp, SNRi_cut=args.snri, flux_lim=args.lim) + Stokes_UV = fits.open(args.file, mode="readonly") + p = pol_map(Stokes_UV, SNRp_cut=args.snrp, SNRi_cut=args.snri, step_vec=args.step_vec, scale_vec=args.scale_vec, flux_lim=args.lim) else: - print("python3 plots.py -f -p -i -l ") + print("python3 plots.py -f -p -i -st -sc -l ")