Tibeuleu/FOC_Reduction
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

Uniform plots and better handled maps for plots

Browse Source
This commit is contained in:
Thibault Barnouin
2024-05-02 12:37:40 +02:00
parent a931e20f5f
commit cbe6aeda45
3 changed files with 73 additions and 63 deletions
Download Patch File Download Diff File Expand all files Collapse all files

12
src/FOC_reduction.py
View File

@@ -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

109
src/lib/plots.py
View File

@@ -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)

15
src/lib/reduction.py
View File

@@ -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)]])