From cbe6aeda4524f93daeb02318d990471df5155f4e Mon Sep 17 00:00:00 2001 From: Thibault Barnouin Date: Thu, 2 May 2024 12:37:40 +0200 Subject: [PATCH] Uniform plots and better handled maps for plots --- src/FOC_reduction.py | 12 ++--- src/lib/plots.py | 109 ++++++++++++++++++++++++------------------- src/lib/reduction.py | 15 +++--- 3 files changed, 73 insertions(+), 63 deletions(-) diff --git a/src/FOC_reduction.py b/src/FOC_reduction.py index 2e60d98..0b029d9 100755 --- a/src/FOC_reduction.py +++ b/src/FOC_reduction.py @@ -18,7 +18,7 @@ from matplotlib.colors import LogNorm def main(target=None, proposal_id=None, infiles=None, output_dir="./data", crop=False, interactive=False): # Reduction parameters # Deconvolution - deconvolve = True + deconvolve = False if deconvolve: # from lib.deconvolve import from_file_psf psf = 'gaussian' # Can be user-defined as well @@ -33,9 +33,9 @@ def main(target=None, proposal_id=None, infiles=None, output_dir="./data", crop= display_crop = False # Background estimation - error_sub_type = 'sturges' # sqrt, sturges, rice, scott, freedman-diaconis (default) or shape (example (51, 51)) - subtract_error = 0.25 - display_bkg = True + error_sub_type = 'freedman-diaconis' # sqrt, sturges, rice, scott, freedman-diaconis (default) or shape (example (51, 51)) + subtract_error = 0.50 + display_bkg = False # Data binning rebin = True @@ -45,12 +45,12 @@ def main(target=None, proposal_id=None, infiles=None, output_dir="./data", crop= # Alignement align_center = 'center' # If None will not align the images - display_align = True + display_align = False display_data = False # Smoothing smoothing_function = 'combine' # gaussian_after, weighted_gaussian_after, gaussian, weighted_gaussian or combine - smoothing_FWHM = 0.10 # If None, no smoothing is done + smoothing_FWHM = 0.20 # If None, no smoothing is done smoothing_scale = 'arcsec' # pixel or arcsec # Rotation diff --git a/src/lib/plots.py b/src/lib/plots.py index fa12277..cca65ee 100755 --- a/src/lib/plots.py +++ b/src/lib/plots.py @@ -90,7 +90,7 @@ def plot_obs(data_array, headers, rectangle=None, savename=None, plots_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']]) shape = np.array((3, nb_obs)) - fig, ax = plt.subplots(shape[0], shape[1], figsize=(10, 10), dpi=200, + 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) @@ -107,7 +107,11 @@ def plot_obs(data_array, headers, rectangle=None, savename=None, plots_folder="" else: ax_curr = ax[r_ax] # plots - vmin, vmax = convert*data[data > 0.].min()/10., convert*data[data > 0.].max() + 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() for key, value in [["cmap", [["cmap", "gray"]]], ["norm", [["norm", LogNorm(vmin, vmax)]]]]: try: _ = kwargs[key] @@ -127,7 +131,7 @@ def plot_obs(data_array, headers, rectangle=None, savename=None, plots_folder="" 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.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}$]") if not (savename is None): @@ -160,6 +164,10 @@ def plot_Stokes(Stokes, savename=None, plots_folder=""): 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 wcs = WCS(Stokes[0]).deepcopy() @@ -244,17 +252,23 @@ def polarisation_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[np.argmax([Stokes[i].header['datatype'] == 'I_stokes' for i in range(len(Stokes))])] - stk_cov = Stokes[np.argmax([Stokes[i].header['datatype'] == 'IQU_cov_matrix' for i in range(len(Stokes))])] - pol = Stokes[np.argmax([Stokes[i].header['datatype'] == 'Pol_deg_debiased' for i in range(len(Stokes))])] - pol_err = Stokes[np.argmax([Stokes[i].header['datatype'] == 'Pol_deg_err' for i in range(len(Stokes))])] - pang = Stokes[np.argmax([Stokes[i].header['datatype'] == 'Pol_ang' for i in range(len(Stokes))])] + 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[np.argmax([Stokes[i].header['datatype'] == 'Data_mask' for i in range(len(Stokes))])].data.astype(bool) + data_mask = Stokes['Data_mask'].data.astype(bool).copy() except KeyError: data_mask = np.ones(stkI.shape).astype(bool) + for dataset in [stkI, pol, pol_err, pang]: + dataset[np.logical_not(data_mask)] = np.nan + for i in range(3): + 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'] wcs = WCS(Stokes[0]).deepcopy() @@ -264,95 +278,94 @@ def polarisation_map(Stokes, data_mask=None, rectangle=None, SNRp_cut=3., SNRi_c plot_Stokes(Stokes, savename=savename, plots_folder=plots_folder) # Compute SNR and apply cuts - poldata, pangdata = pol.data.copy(), pang.data.copy() - maskP = pol_err.data > 0 - SNRp = np.zeros(pol.data.shape) - SNRp[maskP] = pol.data[maskP]/pol_err.data[maskP] + poldata, pangdata = pol.copy(), pang.copy() + maskP = pol_err > 0 + SNRp = np.ones(pol.shape)*np.nan + SNRp[maskP] = pol[maskP]/pol_err[maskP] - maskI = stk_cov.data[0, 0] > 0 - SNRi = np.zeros(stkI.data.shape) - SNRi[maskI] = stkI.data[maskI]/np.sqrt(stk_cov.data[0, 0][maskI]) + maskI = stk_cov[0, 0] > 0 + 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 pangdata[np.logical_not(mask)] = np.nan # Look for pixel of max polarisation - if np.isfinite(pol.data).any(): - p_max = np.max(pol.data[np.isfinite(pol.data)]) - x_max, y_max = np.unravel_index(np.argmax(pol.data == p_max), pol.data.shape) + if np.isfinite(pol).any(): + p_max = np.max(pol[np.isfinite(pol)]) + x_max, y_max = np.unravel_index(np.argmax(pol == p_max), pol.shape) else: print("No pixel with polarisation information above requested SNR.") # Plot the map plt.rcParams.update({'font.size': 10}) plt.rcdefaults() - fig, ax = plt.subplots(figsize=(10, 10), subplot_kw=dict(projection=wcs)) + 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) + # fig.subplots_adjust(hspace=0, wspace=0, left=0.102, right=1.02) if display.lower() in ['intensity']: # If no display selected, show intensity map display = 'i' if flux_lim is None: if mask.sum() > 0.: - vmin, vmax = 1./2.*np.median(np.sqrt(stk_cov.data[0, 0][mask])*convert_flux), np.max(stkI.data[stkI.data > 0.]*convert_flux) + vmin, vmax = 1./2.*np.median(np.sqrt(stk_cov[0, 0][mask])*convert_flux), np.max(stkI[stkI > 0.]*convert_flux) else: - vmin, vmax = 1./2.*np.median(np.sqrt(stk_cov.data[0, 0][stkI.data > 0.])*convert_flux), np.max(stkI.data[stkI.data > 0.]*convert_flux) + vmin, vmax = 1./2.*np.median(np.sqrt(stk_cov[0, 0][stkI > 0.])*convert_flux), np.max(stkI[stkI > 0.]*convert_flux) else: vmin, vmax = flux_lim - im = ax.imshow(stkI.data*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.) 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 print("Total flux contour levels : ", levelsI) - ax.contour(stkI.data*convert_flux, levels=levelsI, colors='grey', linewidths=0.5) + ax.contour(stkI*convert_flux, levels=levelsI, colors='grey', linewidths=0.5) elif display.lower() in ['pol_flux']: # Display polarisation flux display = 'pf' if flux_lim is None: if mask.sum() > 0.: - vmin, vmax = 1./2.*np.median(np.sqrt(stk_cov.data[0, 0][mask])*convert_flux), np.max(stkI.data[stkI.data > 0.]*convert_flux) + vmin, vmax = 1./2.*np.median(np.sqrt(stk_cov[0, 0][mask])*convert_flux), np.max(stkI[stkI > 0.]*convert_flux) else: - vmin, vmax = 1./2.*np.median(np.sqrt(stk_cov.data[0, 0][stkI.data > 0.])*convert_flux), np.max(stkI.data[stkI.data > 0.]*convert_flux) + vmin, vmax = 1./2.*np.median(np.sqrt(stk_cov[0, 0][stkI > 0.])*convert_flux), np.max(stkI[stkI > 0.]*convert_flux) else: vmin, vmax = flux_lim - im = ax.imshow(stkI.data*convert_flux*pol.data, 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.) 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) print("Polarized flux contour levels : ", levelsPf) - ax.contour(stkI.data*convert_flux*pol.data, levels=levelsPf, colors='grey', linewidths=0.5) + ax.contour(stkI*convert_flux*pol, levels=levelsPf, colors='grey', linewidths=0.5) elif display.lower() in ['p', 'pol', 'pol_deg']: # Display polarisation degree map display = 'p' vmin, vmax = 0., 100. - im = ax.imshow(pol.data*100., vmin=vmin, vmax=vmax, aspect='equal', cmap='inferno', alpha=1.) + im = ax.imshow(pol*100., vmin=vmin, vmax=vmax, aspect='equal', cmap='inferno', alpha=1.) fig.colorbar(im, ax=ax, aspect=50, shrink=0.75, pad=0.025, label=r"$P$ [%]") elif display.lower() in ['pa', 'pang', 'pol_ang']: # Display polarisation degree map display = 'pa' vmin, vmax = 0., 180. - im = ax.imshow(princ_angle(pang.data), vmin=vmin, vmax=vmax, aspect='equal', cmap='inferno', alpha=1.) + im = ax.imshow(princ_angle(pang), vmin=vmin, vmax=vmax, aspect='equal', cmap='inferno', alpha=1.) 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']: # Display polarisation degree error map display = 's_p' if (SNRp > SNRp_cut).any(): - vmin, vmax = 0., np.max(pol_err.data[SNRp > SNRp_cut])*100. - p_err = deepcopy(pol_err.data) - p_err[p_err > vmax/100.] = np.nan - im = ax.imshow(p_err*100., vmin=vmin, vmax=vmax, aspect='equal', cmap='inferno', alpha=1.) + 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.) else: - im = ax.imshow(pol_err.data*100., aspect='equal', cmap='inferno', alpha=1.) + vmin, vmax = 0., 100. + im = ax.imshow(pol_err*100., vmin=vmin, vmax=vmax, aspect='equal', cmap='inferno_r', alpha=1.) 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']: # Display intensity error map display = 's_i' if (SNRi > SNRi_cut).any(): - vmin, vmax = 1./2.*np.median(np.sqrt(stk_cov.data[0, 0][stk_cov.data[0, 0] > 0.]) * - convert_flux), np.max(np.sqrt(stk_cov.data[0, 0][stk_cov.data[0, 0] > 0.])*convert_flux) - im = ax.imshow(np.sqrt(stk_cov.data[0, 0])*convert_flux, norm=LogNorm(vmin, vmax), aspect='equal', cmap='inferno', alpha=1.) + 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.) else: - im = ax.imshow(np.sqrt(stk_cov.data[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.) 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']: # Display I_stokes signal-to-noise map @@ -381,15 +394,15 @@ def polarisation_map(Stokes, data_mask=None, rectangle=None, SNRp_cut=3., SNRi_c else: # Defaults to intensity map if mask.sum() > 0.: - vmin, vmax = 1.*np.mean(np.sqrt(stk_cov.data[0, 0][mask])*convert_flux), np.max(stkI.data[stkI.data > 0.]*convert_flux) + vmin, vmax = 1.*np.mean(np.sqrt(stk_cov[0, 0][mask])*convert_flux), np.max(stkI[stkI > 0.]*convert_flux) else: - vmin, vmax = 1.*np.mean(np.sqrt(stk_cov.data[0, 0][stkI.data > 0.])*convert_flux), np.max(stkI.data[stkI.data > 0.]*convert_flux) - im = ax.imshow(stkI.data*convert_flux, norm=LogNorm(vmin, vmax), aspect='equal', cmap='inferno', alpha=1.) + 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.) 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[data_mask].sum() - I_diluted_err = np.sqrt(np.sum(stk_cov.data[0, 0][data_mask])) + 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'] @@ -407,7 +420,7 @@ def polarisation_map(Stokes, data_mask=None, rectangle=None, SNRp_cut=3., SNRi_c poldata[np.isfinite(poldata)] = 1./2. step_vec = 1 vec_scale = 2. - X, Y = np.meshgrid(np.arange(stkI.data.shape[1]), np.arange(stkI.data.shape[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.2, linewidth=0.3, color='w', edgecolor='k') @@ -429,7 +442,7 @@ def polarisation_map(Stokes, data_mask=None, rectangle=None, SNRp_cut=3., SNRi_c # Display instrument FOV if not (rectangle is None): x, y, width, height, angle, color = rectangle - x, y = np.array([x, y]) - np.array(stkI.data.shape)/2. + 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)) @@ -2053,7 +2066,7 @@ class pol_map(object): def submit_save(expression): ax_text_save.set(visible=False) if expression != '': - save_fig, save_ax = plt.subplots(figsize=(12, 10), layout='tight', subplot_kw=dict(projection=self.wcs)) + save_fig, save_ax = plt.subplots(figsize=(12, 10), layout='constrained', subplot_kw=dict(projection=self.wcs)) self.ax_cosmetics(ax=save_ax) self.display(fig=save_fig, ax=save_ax) self.pol_vector(fig=save_fig, ax=save_ax) diff --git a/src/lib/reduction.py b/src/lib/reduction.py index ba7df4a..23551b5 100755 --- a/src/lib/reduction.py +++ b/src/lib/reduction.py @@ -276,7 +276,7 @@ def crop_array(data_array, headers, error_array=None, data_mask=None, step=5, nu if display: plt.rcParams.update({'font.size': 15}) - fig, ax = plt.subplots(figsize=(10, 10)) + fig, ax = plt.subplots(figsize=(10, 10), layout='constrained') convert_flux = headers[0]['photflam'] data = deepcopy(data_array[0]*convert_flux) data[data <= data[data > 0.].min()] = data[data > 0.].min() @@ -301,16 +301,13 @@ def crop_array(data_array, headers, error_array=None, data_mask=None, step=5, nu xycoords='axes fraction') ax.annotate(str(exptime)+" s", color='white', fontsize=10, xy=(0.80, 0.02), xycoords='axes fraction') - ax.set( # title="Location of cropped image.", - xlabel='pixel offset', - ylabel='pixel offset') + ax.set(title="Location of cropped image.", xlabel='pixel offset', ylabel='pixel offset') - fig.subplots_adjust(hspace=0, wspace=0, right=0.85) - cbar_ax = fig.add_axes([0.9, 0.12, 0.02, 0.75]) - fig.colorbar(im, cax=cbar_ax, label=r"Flux [$ergs \cdot cm^{-2} \cdot s^{-1} \cdot \AA^{-1}$]") + # fig.subplots_adjust(hspace=0, wspace=0, right=0.85) + # cbar_ax = fig.add_axes([0.9, 0.12, 0.02, 0.75]) + fig.colorbar(im, ax=ax, label=r"Flux [$ergs \cdot cm^{-2} \cdot s^{-1} \cdot \AA^{-1}$]") if savename is not None: - # fig.suptitle(savename+'_'+filt+'_crop_region') fig.savefig("/".join([plots_folder, savename+'_'+filt+'_crop_region.png']), bbox_inches='tight') plot_obs(data_array, headers, vmin=convert_flux*data_array[data_array > 0.].mean()/5., @@ -1406,7 +1403,7 @@ def rotate_Stokes(I_stokes, Q_stokes, U_stokes, Stokes_cov, data_mask, headers, for i, head in enumerate(headers): ang[i] = -head['orientat'] ang = ang.mean() - alpha = np.radians(ang) + alpha = np.pi/180.*ang mrot = np.array([[1., 0., 0.], [0., np.cos(2.*alpha), np.sin(2.*alpha)], [0, -np.sin(2.*alpha), np.cos(2.*alpha)]])