add the possibility to fix flux limits in output maps
This commit is contained in:
Thibault Barnouin
@@ -63,6 +63,7 @@ def main(target=None, proposal_id=None, infiles=None, output_dir="./data"):
|
||||
# Polarization map output
|
||||
SNRp_cut = 3. #P measurments with SNR>3
|
||||
SNRi_cut = 30. #I measurments with SNR>30, which implies an uncertainty in P of 4.7%.
|
||||
flux_lim = [5e-19,5e-14] #lowest and highest flux displayed on plot, defaults to bkg and maximum in cut if None
|
||||
vec_scale = 2.0
|
||||
step_vec = 1 #plot all vectors in the array. if step_vec = 2, then every other vector will be plotted
|
||||
# if step_vec = 0 then all vectors are displayed at full length
|
||||
@@ -182,19 +183,19 @@ def main(target=None, proposal_id=None, infiles=None, output_dir="./data"):
|
||||
print("PA_bkg = {0:.1f} ± {1:.1f} °".format(PA_bkg[0,0],np.ceil(s_PA_bkg[0,0]*10.)/10.))
|
||||
# Plot polarization map (Background is either total Flux, Polarization degree or Polarization degree error).
|
||||
if px_scale.lower() not in ['full','integrate'] and final_display:
|
||||
proj_plots.polarization_map(deepcopy(Stokes_test), data_mask, SNRp_cut=SNRp_cut, SNRi_cut=SNRi_cut, step_vec=step_vec, vec_scale=vec_scale, savename="_".join([figname,figtype]), plots_folder=plots_folder)
|
||||
proj_plots.polarization_map(deepcopy(Stokes_test), data_mask, SNRp_cut=SNRp_cut, SNRi_cut=SNRi_cut, step_vec=step_vec, vec_scale=vec_scale, savename="_".join([figname,figtype,"I"]), plots_folder=plots_folder, display='Intensity')
|
||||
proj_plots.polarization_map(deepcopy(Stokes_test), data_mask, SNRp_cut=SNRp_cut, SNRi_cut=SNRi_cut, step_vec=step_vec, vec_scale=vec_scale, savename="_".join([figname,figtype,"P_flux"]), plots_folder=plots_folder, display='Pol_Flux')
|
||||
proj_plots.polarization_map(deepcopy(Stokes_test), data_mask, SNRp_cut=SNRp_cut, SNRi_cut=SNRi_cut, step_vec=step_vec, vec_scale=vec_scale, savename="_".join([figname,figtype,"P"]), plots_folder=plots_folder, display='Pol_deg')
|
||||
proj_plots.polarization_map(deepcopy(Stokes_test), data_mask, SNRp_cut=SNRp_cut, SNRi_cut=SNRi_cut, step_vec=step_vec, vec_scale=vec_scale, savename="_".join([figname,figtype,"PA"]), plots_folder=plots_folder, display='Pol_ang')
|
||||
proj_plots.polarization_map(deepcopy(Stokes_test), data_mask, SNRp_cut=SNRp_cut, SNRi_cut=SNRi_cut, step_vec=step_vec, vec_scale=vec_scale, savename="_".join([figname,figtype,"I_err"]), plots_folder=plots_folder, display='I_err')
|
||||
proj_plots.polarization_map(deepcopy(Stokes_test), data_mask, SNRp_cut=SNRp_cut, SNRi_cut=SNRi_cut, step_vec=step_vec, vec_scale=vec_scale, savename="_".join([figname,figtype,"P_err"]), plots_folder=plots_folder, display='Pol_deg_err')
|
||||
proj_plots.polarization_map(deepcopy(Stokes_test), data_mask, SNRp_cut=SNRp_cut, SNRi_cut=SNRi_cut, step_vec=step_vec, vec_scale=vec_scale, savename="_".join([figname,figtype,"SNRi"]), plots_folder=plots_folder, display='SNRi')
|
||||
proj_plots.polarization_map(deepcopy(Stokes_test), data_mask, SNRp_cut=SNRp_cut, SNRi_cut=SNRi_cut, step_vec=step_vec, vec_scale=vec_scale, savename="_".join([figname,figtype,"SNRp"]), plots_folder=plots_folder, display='SNRp')
|
||||
proj_plots.polarization_map(deepcopy(Stokes_test), data_mask, SNRp_cut=SNRp_cut, SNRi_cut=SNRi_cut, flux_lim=flux_lim, step_vec=step_vec, vec_scale=vec_scale, savename="_".join([figname,figtype]), plots_folder=plots_folder)
|
||||
proj_plots.polarization_map(deepcopy(Stokes_test), data_mask, SNRp_cut=SNRp_cut, SNRi_cut=SNRi_cut, flux_lim=flux_lim, step_vec=step_vec, vec_scale=vec_scale, savename="_".join([figname,figtype,"I"]), plots_folder=plots_folder, display='Intensity')
|
||||
proj_plots.polarization_map(deepcopy(Stokes_test), data_mask, SNRp_cut=SNRp_cut, SNRi_cut=SNRi_cut, flux_lim=flux_lim, step_vec=step_vec, vec_scale=vec_scale, savename="_".join([figname,figtype,"P_flux"]), plots_folder=plots_folder, display='Pol_Flux')
|
||||
proj_plots.polarization_map(deepcopy(Stokes_test), data_mask, SNRp_cut=SNRp_cut, SNRi_cut=SNRi_cut, flux_lim=flux_lim, step_vec=step_vec, vec_scale=vec_scale, savename="_".join([figname,figtype,"P"]), plots_folder=plots_folder, display='Pol_deg')
|
||||
proj_plots.polarization_map(deepcopy(Stokes_test), data_mask, SNRp_cut=SNRp_cut, SNRi_cut=SNRi_cut, flux_lim=flux_lim, step_vec=step_vec, vec_scale=vec_scale, savename="_".join([figname,figtype,"PA"]), plots_folder=plots_folder, display='Pol_ang')
|
||||
proj_plots.polarization_map(deepcopy(Stokes_test), data_mask, SNRp_cut=SNRp_cut, SNRi_cut=SNRi_cut, flux_lim=flux_lim, step_vec=step_vec, vec_scale=vec_scale, savename="_".join([figname,figtype,"I_err"]), plots_folder=plots_folder, display='I_err')
|
||||
proj_plots.polarization_map(deepcopy(Stokes_test), data_mask, SNRp_cut=SNRp_cut, SNRi_cut=SNRi_cut, flux_lim=flux_lim, step_vec=step_vec, vec_scale=vec_scale, savename="_".join([figname,figtype,"P_err"]), plots_folder=plots_folder, display='Pol_deg_err')
|
||||
proj_plots.polarization_map(deepcopy(Stokes_test), data_mask, SNRp_cut=SNRp_cut, SNRi_cut=SNRi_cut, flux_lim=flux_lim, step_vec=step_vec, vec_scale=vec_scale, savename="_".join([figname,figtype,"SNRi"]), plots_folder=plots_folder, display='SNRi')
|
||||
proj_plots.polarization_map(deepcopy(Stokes_test), data_mask, SNRp_cut=SNRp_cut, SNRi_cut=SNRi_cut, flux_lim=flux_lim, step_vec=step_vec, vec_scale=vec_scale, savename="_".join([figname,figtype,"SNRp"]), plots_folder=plots_folder, display='SNRp')
|
||||
elif final_display:
|
||||
proj_plots.polarization_map(deepcopy(Stokes_test), data_mask, SNRp_cut=SNRp_cut, SNRi_cut=SNRi_cut, savename="_".join([figname,figtype]), plots_folder=plots_folder, display='integrate')
|
||||
elif px_scale.lower() not in ['full', 'integrate']:
|
||||
pol_map = proj_plots.pol_map(Stokes_test, SNRp_cut=SNRp_cut, SNRi_cut=SNRi_cut)
|
||||
pol_map = proj_plots.pol_map(Stokes_test, SNRp_cut=SNRp_cut, SNRi_cut=SNRi_cut, flux_lim=flux_lim)
|
||||
|
||||
return 0
|
||||
|
||||
|
||||
@@ -3,11 +3,12 @@ from getopt import getopt, error as get_error
|
||||
from sys import argv
|
||||
|
||||
arglist = argv[1:]
|
||||
options = "hf:p:i:"
|
||||
long_options = ["help","fits=","snrp=","snri="]
|
||||
options = "hf:p:i:l:"
|
||||
long_options = ["help","fits=","snrp=","snri=","lim="]
|
||||
|
||||
fits_path = None
|
||||
SNRp_cut, SNRi_cut = 3, 30
|
||||
flux_lim = None
|
||||
out_txt = None
|
||||
|
||||
try:
|
||||
@@ -15,13 +16,15 @@ try:
|
||||
|
||||
for curr_arg, curr_val in arg:
|
||||
if curr_arg in ("-h", "--help"):
|
||||
print("python3 analysis.py -f <path_to_reduced_fits> -p <SNRp_cut> -i <SNRi_cut>")
|
||||
print("python3 analysis.py -f <path_to_reduced_fits> -p <SNRp_cut> -i <SNRi_cut> -l <flux_lim>")
|
||||
elif curr_arg in ("-f", "--fits"):
|
||||
fits_path = str(curr_val)
|
||||
elif curr_arg in ("-p", "--snrp"):
|
||||
SNRp_cut = int(curr_val)
|
||||
elif curr_arg in ("-i", "--snri"):
|
||||
SNRi_cut = int(curr_val)
|
||||
elif curr_arg in ("-l", "--lim"):
|
||||
flux_lim = list("".join(curr_val).split(','))
|
||||
except get_error as err:
|
||||
print(str(err))
|
||||
|
||||
@@ -30,7 +33,7 @@ if not fits_path is None:
|
||||
from lib.plots import pol_map
|
||||
|
||||
Stokes_UV = fits.open(fits_path)
|
||||
p = pol_map(Stokes_UV, SNRp_cut=SNRp_cut, SNRi_cut=SNRi_cut)
|
||||
p = pol_map(Stokes_UV, SNRp_cut=SNRp_cut, SNRi_cut=SNRi_cut,flux_lim=flux_lim)
|
||||
|
||||
else:
|
||||
print("python3 analysis.py -f <path_to_reduced_fits> -p <SNRp_cut> -i <SNRi_cut>")
|
||||
print("python3 analysis.py -f <path_to_reduced_fits> -p <SNRp_cut> -i <SNRi_cut> -l <flux_lim>")
|
||||
|
||||
@@ -219,7 +219,7 @@ def plot_Stokes(Stokes, savename=None, plots_folder=""):
|
||||
|
||||
|
||||
def polarization_map(Stokes, data_mask=None, rectangle=None, SNRp_cut=3., SNRi_cut=30.,
|
||||
step_vec=1, vec_scale=2., savename=None, plots_folder="", display="default"):
|
||||
flux_lim=None, step_vec=1, vec_scale=2., savename=None, plots_folder="", display="default"):
|
||||
"""
|
||||
Plots polarization map from Stokes HDUList.
|
||||
----------
|
||||
@@ -239,6 +239,10 @@ def polarization_map(Stokes, data_mask=None, rectangle=None, SNRp_cut=3., SNRi_c
|
||||
Cut that should be applied to the signal-to-noise ratio on I.
|
||||
Any SNR < SNRi_cut won't be displayed.
|
||||
Defaults to 30. This value implies an uncertainty in P of 4.7%
|
||||
flux_lim : float list, optional
|
||||
Limits that should be applied to the flux colorbar.
|
||||
Defaults to None, limits are computed on the background value and the
|
||||
maximum value in the cut.
|
||||
step_vec : int, optional
|
||||
Number of steps between each displayed polarization vector.
|
||||
If step_vec = 2, every other vector will be displayed.
|
||||
@@ -320,10 +324,13 @@ def polarization_map(Stokes, data_mask=None, rectangle=None, SNRp_cut=3., SNRi_c
|
||||
if display.lower() in ['intensity']:
|
||||
# If no display selected, show intensity map
|
||||
display='i'
|
||||
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)
|
||||
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)
|
||||
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)
|
||||
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 = flux_lim
|
||||
im = ax.imshow(stkI.data*convert_flux, norm=LogNorm(vmin,vmax), aspect='equal', cmap='inferno', alpha=1.)
|
||||
cbar = plt.colorbar(im, cax=cbar_ax, label=r"$F_{\lambda}$ [$ergs \cdot cm^{-2} \cdot s^{-1} \cdot \AA^{-1}$]")
|
||||
levelsI = np.linspace(vmax*0.01, vmax*0.99, 10)
|
||||
@@ -334,10 +341,13 @@ def polarization_map(Stokes, data_mask=None, rectangle=None, SNRp_cut=3., SNRi_c
|
||||
# Display polarisation flux
|
||||
display='pf'
|
||||
pf_mask = (stkI.data > 0.) * (pol.data > 0.)
|
||||
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)
|
||||
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)
|
||||
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)
|
||||
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 = flux_lim
|
||||
im = ax.imshow(stkI.data*convert_flux*pol.data, norm=LogNorm(vmin,vmax), aspect='equal', cmap='inferno', alpha=1.)
|
||||
cbar = plt.colorbar(im, cax=cbar_ax, 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)
|
||||
@@ -1353,13 +1363,14 @@ class pol_map(object):
|
||||
"""
|
||||
Class to interactively study polarization maps.
|
||||
"""
|
||||
def __init__(self, Stokes, SNRp_cut=3., SNRi_cut=30., selection=None):
|
||||
def __init__(self, Stokes, SNRp_cut=3., SNRi_cut=30., flux_lim=None, selection=None):
|
||||
|
||||
if type(Stokes) == str:
|
||||
Stokes = fits.open(Stokes)
|
||||
self.Stokes = deepcopy(Stokes)
|
||||
self.SNRp_cut = SNRp_cut
|
||||
self.SNRi_cut = SNRi_cut
|
||||
self.flux_lim = flux_lim
|
||||
self.SNRi = deepcopy(self.SNRi_cut)
|
||||
self.SNRp = deepcopy(self.SNRp_cut)
|
||||
self.region = None
|
||||
@@ -1782,18 +1793,26 @@ class pol_map(object):
|
||||
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})
|
||||
ax.add_artist(self.north_dir)
|
||||
|
||||
def display(self, fig=None, ax=None):
|
||||
def display(self, fig=None, ax=None, flux_lim=None):
|
||||
norm = None
|
||||
if self.display_selection is None:
|
||||
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.convert_flux
|
||||
vmin, vmax = 1./2.*np.median(self.data[self.data > 0.]), np.max(self.data[self.data > 0.])
|
||||
if flux_lim is None:
|
||||
vmin, vmax = 1./2.*np.median(self.data[self.data > 0.]), np.max(self.data[self.data > 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.convert_flux*self.P
|
||||
vmin, vmax = 1./2.*np.median(self.I[self.I > 0.]*self.convert_flux), np.max(self.I[self.I > 0.]*self.convert_flux)
|
||||
if flux_lim is None:
|
||||
vmin, vmax = 1./2.*np.median(self.I[self.I > 0.]*self.convert_flux), np.max(self.I[self.I > 0.]*self.convert_flux)
|
||||
else:
|
||||
vmin, vmax = flux_lim
|
||||
norm = LogNorm(vmin, vmax)
|
||||
label = r"$F_{\lambda} \cdot P$ [$ergs \cdot cm^{-2} \cdot s^{-1} \cdot \AA^{-1}$]"
|
||||
elif self.display_selection.lower() in ['pol_deg']:
|
||||
@@ -1853,11 +1872,11 @@ class pol_map(object):
|
||||
ax = self.ax
|
||||
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.1, linewidth=0.5, color='white',edgecolor='black')
|
||||
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.15, linewidth=0.5, 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.1, linewidth=0.5, color='white',edgecolor='black')
|
||||
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.15, linewidth=0.5, color='white',edgecolor='black')
|
||||
fig.canvas.draw_idle()
|
||||
|
||||
def pol_int(self, fig=None, ax=None):
|
||||
|
||||
Reference in New Issue
Block a user